Extremal correlation coefficient for functional data

Mihyun Kim$^1$ and  Piotr Kokoszka$^2$

${}^1$Department of Statistics, West Virginia University, Morgantown, WV, USA

${}^2$Department of Statistics, Colorado State University,
Fort Collins, CO, USA

Address for correspondence: Mihyun Kim,
Department of Statistics,
West Virginia University, Morgantown, WV 26506, USA.
Email<EMAIL_ADDRESS>
We propose a coefficient that measures dependence in
paired samples of functions. It has properties similar to the
Pearson correlation,  but  differs in significant ways:
1) it is designed to measure dependence between curves, 2)
it focuses only on extreme curves. The new coefficient is derived
within the framework of regular variation in Banach
spaces. A consistent estimator is proposed and justified
by an asymptotic analysis and a simulation study.
The usefulness of the new coefficient is illustrated on
financial and and climate functional data.

Keywords: correlation, extremes, functional data.

§ INTRODUCTION

Due to the growing  impact of extreme events related, for example,
to financial downturns or unusual weather, there has been increasing
interest in developing statistical tools to study patterns of extreme curves.
This is to a large extent due to the increasing availability of high resolution
data; time series of asset prices can be constructed at practically
any temporal resolution, and modern weather databases contain
measurements  at hourly or even higher frequencies. Such data
can be interpreted as curves, e.g., one curve per day,
and provide  more complete information than a single number
per day, like the total return or the maximum temperature.
We propose a coefficient that quantifies the tendency of paired extreme curves to exhibit similar patterns simultaneously. Two examples of the type of questions that the tool deals with are the following: 1) During a stock market crisis, such as the market decline  due to the COVID-19 pandemic, do
returns  of  different sectors of the economy exhibit similar extreme daily trajectories? 2) How likely is location A to experience  a similar
daily pattern of temperature  as location B (on the same day)  during a heat
wave? The coefficient we propose  focuses not just on extreme total
return or extreme maximum daily temperature,
but also on the shapes of extreme curves.

There has been some research focusing on probabilistic and
statistical methods for extreme curves.
Extreme value theory in the  space of continuous functions is
studied in Chapters 9 and 10 of  [6].
Principal component analysis of extreme curves has been studied by [20], [18], and [2]. Extremal properties of scores of functional data were studied by [21] and Kim and Kokoszka ([16], [17]).
Additional,  more closely related papers are introduced as we
develop our approach.
We propose a method for quantifying extremal dependence of paired functional samples, for which there are currently no appropriate tools.

We note that there has been considerable research aimed at quantifying extremal dependence for heavy-tailed random vectors. Ledford and Tawn ([23], [24],
[25]) introduced the coefficient of tail dependence, which was later generalized to the extremogram by
[5]. The extremal dependence measure based on the angular measure of a regularly varying random vector was introduced by [29] and further investigated by
[22]. [15] recently introduced a unified approach for representing tail dependence using random exceedence sets. Those measures for extremes are  designed for random vectors in a Euclidean space. Thus applying any such measures to functional data requires some sort of dimension reduction, e.g., principal component analysis, or data compression like converting daily temperature curves to daily average or maximum values. The reduced data are then analyzed using those tools for multivariate extremes, see, e.g., [27], [7], and [17]. This approach is convenient, but it does not fully utilize all relevant information that functional data contain.

We develop a new measure, the extremal correlation coefficient, that captures the extremal dependence of paired functional samples utilizing the information in the sizes and shapes of the curves. The measure is constructed by computing an inner product of pairs of extreme curves. This approach is closely related to the concept of cosine similarity that is often used, for example, to quantify document similarity in text analysis. While the cosine similarity computes the inner product between two vectors to see whether they are pointing in the same direction using all available pairs,  our measure calculates the inner product between paired extreme curves to see whether they look alike in extremes. Similar ideas have been applied in non-extreme contexts of functional data analysis. [9] introduced a measure, called dynamical correlation, that computes the inner product of all pairs of standardized curves. The concept was further studied by
[31] where an  autocorrelation measure, termed spherical autocorrelation, for functional time series was proposed. These measures are however computed based on the total body of functional data, and so are not suitable for describing extremal dependence.

The coefficient  we develop quantifies extremal dependence in pairs of heavy-tailed functional observations. It is conceptually appealing, as it shares desirable features with the classic correlation coefficient: 1) its values range from -1 to 1, 2) it measures the strength and direction of linear relationship between two extreme curves, 3) if the extremal behavior of two curves is independent, the coefficient is zero. Moreover, it can be used in practice with a relatively simple numerical implementation. We thus hope that such interpretable and tractable tool makes a useful contribution. No measures of extremal dependence for pairs of curves are currently available.

Turning to mathematical challenges, the concept of vague convergence, see e.g., Chapters 2 and 6 of [30], cannot be readily used. The vague convergence, which now provides a standard mathematical framework for extremes in Euclidean spaces, can be defined only on locally compact spaces. Since every locally compact Banach space has finite dimension, a different framework must be used for functional data in Hilbert spaces. We use the theory of regularly varying measures developed by [14] who introduced the notion of $M_0$ convergence, which works for regularly varying measures on complete separable metric spaces. The $M_0$ convergence is further studied by [27], where it is applied to regularly varying time series in a separable Banach space. Within this framework,
to establish the  the consistency of the estimator we propose.
To do so, we proceed through a
a number of $M_0$ convergence results that allow us to apply an abstract Bernstein-type inequality. A method of computing the extremal
correlation coefficients analytically (in relatively simple cases)
is also developed.

The remainder of the paper is organized as follows. In Section <ref>, we review regularly varying random elements in Banach spaces. In Section <ref>,  we extend the concept of regular variation to bivariate random elements in a product Banach space. The extremal correlation coefficient is introduced in Section <ref>, where its asymptotic properties are also studied. Section <ref> contains  a simulation study, and Section <ref> illustrates applications to  intraday return curves  and daily temperature curves.

Theoretical justification of  our approach requires more detailed background and some technical derivations, which are placed in the
online Supplementary Material. Preliminary results for the proof of Theorem <ref> are presented in Section <ref>, followed by its proof  in Section <ref>. In Section <ref>, the proof of Lemma <ref> is presented.

§ REGULAR VARIATION IN BANACH SPACES

This section presents background needed to understand the development in Sections <ref> and <ref>. In Functional Data Analysis, observations are typically treated as elements of $L^2 := L^2(\cT)$, where the measure space $\cT$ is such that $L^2(\cT)$ is a separable Hilbert
space, equipped with the usual inner product $\lip x,y\rip = \int_{\cT} x(t)y(t) dt$.  The $L^2$-norm is then $\| x\| = \lip x, x\rip^{1/2}$ = $ \lp \int_{\cT} x(t)^2 dt \rp^{1/2}$. An introduction to Functional Data Analysis is presented in [19], while a detailed mathematical treatment is given in [13]. While we refer to the elements of
$L^2$ as curves, due to the  examples we consider, the set
$\cT$ can  be a fairly abstract space (a metric Polish space),
for example a spatial domain.

An extreme curve in $L^2$ is defined as a functional object whose size, measured by the $L^2$-norm, is large. The norm can be large for various reasons as long as the area under the squared curve on $\cT$ is large. For example, curves that are far away from the sample mean or that fluctuate a lot around the sample mean will be extreme according to this definition. Extreme functional observations are thus very different from extreme scalar or multivariate observations because there is a multitude of ways in which a curve can be extreme. We informally call functional data heavy-tailed if the probability that an  extreme curve occurs is relatively large.

We now briefly review the $M_0$ convergence in a separable Banach space $\mbB$. In what follows, $\bzero$ is the zero element. Fix a norm $\| \cdot \|_{\mbB}$ and let $B_{\eg} := \{z \in \mbB: \| z\|_{\mbB} <\eg\}$ be the open ball of radius $\eg>0$ centered at the origin. A Borel measure $\mu$ defined on $\mbB_0:=\mbB  \setminus \{\bzero\}$ is said to be boundedly finite if $\mu(A)<\infty$, for all Borel sets that are bounded away from $\bzero$, i.e., $A \cap B_{\eg} = \emptyset$, for some $\eg>0$. Let $M_0 (\mbB)$ be the collection of all such measures on $\mbB_0$. For $\mu_n, \mu \in M_0(\mbB)$, the sequence of $\mu_n$ converges to $\mu$ in the $M_0$ topology ($\mu_n \stackrel{M_0}{\longrightarrow} \mu$),
if $\mu_n(A) \to \mu(A)$, for all bounded away from $\bzero$, $\mu$–continuity Borel sets $A$, i.e., those with $\mu(\partial A)=0$, where $\partial A$ is the boundary of $A$. Equivalently, $\mu_n \stackrel{M_0}{\longrightarrow} \mu$, if $\int_{\mbB} f(x) \mu_n(dx) \to \int_{\mbB} f(x) \mu(dx)$ for all $f \in \cC_0(\mbB)$, where $\cC_0 (\mbB)$ is the class of bounded and continuous functions $f:\mbB_0 \to \mbR$ that vanish on a neighborhood of $\bzero$.

We now define regular variation for random elements in $\mbB$, see Theorem 3.1 of [14] and Chapter 2 of [27]. This concept formalizes the idea of heavy-tailed observations in infinite dimensional spaces.

A random element $X$ in $\mbB$ is regularly varying with index $-\ag$, $\ag> 0$, if there exist a sequence $b(n) \to \infty$ and a measure $\mu$ in $M_0(\mbB)$ such that

\begin{equation} \label{e:def1-H}
nP \lp  \frac{X}{b(n)}  \in \cdot \rp
\stackrel{M_0}{\longrightarrow} \mu, \ \ \ \ n \to \infty,
\end{equation}

where the exponent measure $\mu$ satisfies $\mu(tA) = t^{-\ag}\mu(A)$ for Borel sets $A \subset \mbB_0$.

A possible choice for $b(n)$ is the quantile function, defined by $P(\| X  \|_{\mbB} >b(n)) = n^{-1}$. Roughly speaking, the tail probability of $X$ decays like a power function, $P(\| X \|_{\mbB} > t ) \approx C t^{-\ag}$, as $t\to \infty$.
The following lemma, see [14],
states an equivalent definition of a regularly varying  element in $\mbB$.

A random element $X$ in $\mbB$ is regularly varying with index $-\ag$, $\ag> 0$, if and only if there exist a sequence $b^{\prime}(n) \to \infty$ and a probability measure $\Gg$ on $\mbS:=\{x \in \mbB : \|x \|_{\mbB}=1\}$  (called the angular measure) such that for any $y>0$,

\begin{equation} \label{e:def2-H}
nP \lp \| X  \|_{\mbB} >b^{\prime}(n)y, X/\| X  \|_{\mbB} \in \cdot \rp
\stackrel{w}{\longrightarrow} c y^{-\ag}\Gg, \ \ \ \ n \to \infty,
\end{equation}

for some $c>0$.

The first three orthonormal basis elements in $L^2[0,1]$ defined in basis (left-most); simulated data when $\Gg$ concentrates on $\phi_1$(second from the left); on $\phi_2$ (third left); on $\phi_3$ (fourth left).

If Definition <ref> (or condition def2-H) holds,
we write $X \in RV(-\ag, \Gg)$. The polar representation def2-H provides an intuitive interpretation of regular variation in $\mbB$. It characterizes regular variation of $X$ in $\mbB$ using two components, the tail index $\ag$ and the angular probability measure $\Gg$. The tail index $\ag$ quantifies how heavy the tail distribution of $\|X\|_{\mbB}$  is, e.g., the probability of extreme curves occurring gets higher as $\ag$ gets smaller. While the tail index $\ag$ determines the frequency of occurrence of extreme curves, the angular measure $\Gg$, defined on the unit sphere $\mbS$, fully characterizes the distribution of the shape of the scaled extreme curves, $X/\|X\|_{\mbB}$. To illustrate this, consider a set of orthonormal functions in $L^2([0,1])$ of the form

\begin{equation} \label{e:basis}
\phi_j(t) = \sqrt{2} \sin \lp \lp j-\frac{1}{2}\rp \pi t\rp, \ \ \ \ j=1,2,\ldots, \ \ t \in [0,1].
\end{equation}

The first three functions are shown in the left-most plot of Figure <ref>. We consider a finite-dimensional subspace of $L^2([0,1])$, spanned by the first 9 $\phi_j$'s, for the purpose of simulations. The data generating process is
$X(t) = \sum_{j=1}^{9} Z_{j} \phi_j(t)$,
where $\bZ=[Z_1, \ldots, Z_9]$ is a 9-dimensional random vector with independent components. Suppose that $Z$ is a random variable following a Pareto distribution with tail index $\ag=3$ and $N$ is a normal random variable with mean 0 and variance .5. We consider the following three cases for $\bZ$:

1. $\bZ = [Z, N, N, N, \ldots, N]^{\top}$; the angular measure $\Gg$ concentrates on $\phi_1$.

2. $\bZ = [N, Z, N, N \ldots, N]^{\top}$; the angular measure $\Gg$ concentrates on $\phi_2$.

3. $\bZ = [N, N, Z, N \ldots, N]^{\top}$; the angular measure $\Gg$ concentrates on $\phi_3$.

Figure <ref> displays simulated data with sample size of 100 for each of the three cases. The plots of simulated data clearly show that the angular measure $\Gg$ represents the distribution of the shapes of extreme curves in that they are dominated by the shape of the functional axis $\phi_j$ on which $\Gg$ concentrates.

§ BIVARIATE REGULAR VARIATION IN BANACH SPACES

In order to describe the extremal dependence of two regularly varying random elements $X$ and $Y$ in $L^2$, we need to identify their joint probabilistic behavior. We again study it in the more general space $\mbB^2$. We propose the following definition.

A bivariate random element $[X,Y]^{\top
}$ in $\mbB^2$ is said to be jointly regularly varying with index $-\ag$, $\ag> 0$, if there exist a sequence $b(n) \to \infty$ and a measure $\mu$ in $M_0(\mbB^2)$ such that

\begin{equation} \label{e:def-RVXY}
nP \lp \frac{(X,Y)}{b(n)}  \in \cdot \rp
\stackrel{M_0}{\longrightarrow} \mu, \ \ \ \ n \to \infty,
\end{equation}

where the joint exponent measure $\mu$ satisfies $\mu(tA) = t^{-\ag}\mu(A)$ for Borel sets $A \subset \mbB_0^2$.

We assume that one-dimensional marginal distributions of $\mu$ are non-degenerate, i.e., $\mu_X:=\mu(\cdot \times \mbB)$ and $\mu_Y:=\mu(\mbB \times \cdot)$ are measures in $M_0(\mbB)$ satisfying analogs of def1-H. Since $X$ and $Y$ are normalized by the same function $b(n)$, the marginal distributions are tail equivalent. A possible choice for $b(n)$ is the quantile function, defined by

\[
P(\|(X,Y)\|_{\mbB^2} >b(n)) = n^{-1}.
\]

With this choice, we have that

\[
nP \lp \frac{(X,Y)}{b(n)}  \in \cA_1 \rp = \frac{P \lp (X,Y) \in b(n)\cA_1   \rp }{P(\|(X,Y)\|_{\mbB^2} >b(n)) } =\frac{P \lp (X,Y) \in \cA_{b(n)}  \rp }{P \lp (X,Y) \in \cA_{b(n)} \rp}=1, \]

where $\cA_{r}$ is defined by

\begin{equation} \label{e:A}
\mathcal A_r =\{(x,y)\in \mbB^2: \|(X,Y)\|_{\mbB^2} \ge r\}, \ \ \ \ r>0.
\end{equation}

Thus, it follows from the $M_0$ convergence in def-RVXY and Lemma <ref> that

\begin{equation} \label{e:mu1}
\mu(\cA_1 )=\mu\{(x,y) \in \mbB^2: \|(X,Y)\|_{\mbB^2} >1\}=1,
\end{equation}

which implies that $\mu$ is a probability measure on $\cA_1$. Throughout the paper, we set

\[
\|(x,y)\|_{\mbB^2} := \|x\|_{\mbB} \vee \|y\|_{\mbB}.
\]

This choice of norm works well with the extremal correlation coefficient defined in Section <ref>.

In order to derive the joint angular probability measure of $X$ and $Y$, we consider the polar coordinate transformation $T:\mbB_0^2 \to \lp [0,\infty)^2  \setminus  \{\bzero\} \rp \times \mbS^2$, defined by

\begin{equation}\label{e:T}
T(x,y) = \lp \|x\|_{\mbB}, \|y\|_{\mbB}, \frac{x}{\|x\|_{\mbB}}, \frac{y}{\|y\|_{\mbB}} \rp=: (r_X, r_Y, \thg_X, \thg_Y), \ \ \ \  (x,y) \in \mbB_0^2.
\end{equation}

Using $T$, we obtain an equivalent formulation for a regularly varying random element in $\mbB^2$.

A bivariate random element $[X,Y]^{\top}$ in $\mbB^2$ is regularly varying with index $-\ag$, $\ag>0$, if and only if there exists an exponent measure $\nu$ in $M_0([0,\infty)^2\setminus\{\bzero\})$ and an angular probability measure $\Gg$ in $M_0(\mbS^2)$ such that

\begin{equation} \label{e:pol-RVXY}
nP \lp \frac{(\|X\|_{\mbB}, \|Y\|_{\mbB})}{b(n)}  \in \cdot, \lp X/\|X\|_{\mbB}, Y/\|Y\|_{\mbB} \rp\in \cdot \rp
\stackrel{M_0}{\longrightarrow} \nu \times \Gg, \ \ n \to \infty,
\end{equation}

where $b(n)$ is the increasing sequence in def-RVXY. (We note that $\mu$ in def-RVXY is a measure on $\mbB_0^2$, while $\nu$ in pol-RVXY is a measure on $[0,\infty)^2 \setminus \{\bzero\}$.)

We only show that def-RVXY implies pol-RVXY  since showing the converse is similar. Take any $f \in \cC_0([0,\infty)^2 \times \mbS^2)$. It then follows from the change of variables that

\begin{align*}
& \int_{[0,\infty)^2} \int_{\mbS^2} f(r_X, r_Y, \thg_X, \thg_Y) \  nP \lp \frac{(\|X\|_{\mbB}, \|Y\|_{\mbB})}{b(n)} \in (dr_X, dr_Y), \lp \frac{X}{\|X\|_{\mbB}}, \frac{Y}{\|Y\|_{\mbB}} \rp  \in (d\thg_X, d\thg_Y)\rp  \\
% & = \int_{[0,\infty)^2 \times \mbS^2} f(r_X, r_Y, \thg_X, \thg_Y) \ \frac{1}{V(u)}P \lp \lp u^{-1}X, u^{-1}Y \rp \rp \circ T^{-1} (dr_X, dr_Y, d\thg_X, d\thg_Y) \\
&= \int_{\mbB^2} f(T(x,y)) \ nP \lp \frac{X}{b(n)} \in dx, \frac{Y}{b(n)} \in dy\rp.
\end{align*}

Since $f \in \cC_0([0,\infty)^2 \times \mbS^2)$, there exists a set $A$, bounded away from $\bzero$ in $[0,\infty)^2 \times \mbS^2$, such that $f(T(x,y)) =0$ if $T(x,y) \notin A$. Then we have that $f(T(x,y)) =0$ if $(x,y) \notin T^{-1}(A)$. Since $T^{-1}(A)$ is bounded away from $\bzero$ in $\mbB^2$, we have that $f \circ T \in \cC_0(\mbB^2)$. Then by the $M_0$ convergence in Definition <ref>, we have that

\begin{align*}
&\int_{\mbB^2} (f \circ T)(x,y) \ nP \lp \frac{X}{b(n)} \in dx, \frac{Y}{b(n)} \in dy\rp \\
&\to \int_{\mbB^2} (f \circ T)(x,y) \mu(dx,dy) = \int_{T(\mbB^2)} f(r_X, r_Y, \thg_X, \thg_Y) \mu\circ T^{-1}(dr_X, dr_Y, d\thg_X, d\thg_Y).
% &=c\int_{\mbS^2} \int_{[0,\infty)^2} f(r_X, r_Y, \thg_X, \thg_Y) \nu(dr_X, dr_Y) \Gg (d\thg_X, d\thg_Y)
\end{align*}

To investigate the form of $\mu\circ T^{-1}$, take any $t>0$ and Borel set $S \subset \mbS^2$. It then follows from the homogeneity property of $\mu$ that

\begin{align*}
&\mu\circ T^{-1} ( [0,\infty)^2\setminus[0,t]^2 \times  S) \\
&=  \mu \lbr  (x,y) \in \mbB^2_0: \|x\|_{\mbB} \vee  \|y\|_{\mbB} >t, (x/\|x\|_{\mbB}, y/\|y\|_{\mbB}) \in S  \rbr  \\
&=  \mu \lbr  (x,y) \in \mbB^2_0: \|x\|_{\mbB} \vee  \|y\|_{\mbB} >t  \rbr \times  \\
&\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \   \frac{ t^{-\ag}\mu \lbr  (x,y) \in \mbB^2_0: \|x\|_{\mbB} \vee  \|y\|_{\mbB} >1, (x/\|x\|_{\mbB}, y/\|y\|_{\mbB}) \in S  \rbr}{  t^{-\ag}\mu \lbr  (x,y) \in \mbB^2_0: \|x\|_{\mbB} \vee  \|y\|_{\mbB} >1 \rbr}.
\end{align*}

It then follows from mu1 that

\[
\mu\circ T^{-1} ( [0,\infty)^2\setminus[0,t]^2 \times  S) = \nu([0,\infty)^2\setminus[0,t]^2) \Gg(S),
\]

\begin{align*}
&\nu(A) := \mu \lbr (x,y) \in \mbB^2_0: (\|x\|_{\mbB},  \|y\|_{\mbB}) \in A \rbr,\ \ \ \ A \subset [0,\infty)^2 \setminus \{\bzero\};\\
&\Gg(S) := \mu \lbr (x,y) \in \mbB^2_0: \|x\|_{\mbB} \vee  \|y\|_{\mbB} > 1, (x/\|x\|_{\mbB}, y/\|y\|_{\mbB}) \in S  \rbr, \ \ \ \ S\subset \mbS^2.
\end{align*}

Thus, $\mu \circ T^{-1}$ has the product form such that on $([0,\infty)^2\setminus \{\bzero\}) \times \mbS^2$

\begin{equation} \label{e:nuGg}
\mu \circ T^{-1} = \nu\times \Gg,
\end{equation}

which completes the proof.

Convergence pol-RVXY can be understood as a polar representation of the bivariate regular variation of $[X,Y]^{\top}$ in $\mbB^2$.
The difference between $\nu$ in pol-RVXY and $\mu$ in def-RVXY is that $\mu$ describes the joint behavior of extreme curves $X$ and $Y$ in $\mbB_0^2$, but $\nu$ describes extremal dependence between the sizes $\|X\|_{\mbB}$ and $\|Y\|_{\mbB}$ in $[0,\infty)^2\setminus \bzero$.
The joint behavior of $X$ and $Y$ in extremes is thus characterized by two measures, $\nu$ on $[0,\infty)^2  \setminus  \{\bzero\}$  and $\Gg$ on $\mbS^2$. The measure $\nu$ describes the joint behavior of $\|X\|_{\mbB}$ and $\|Y\|_{\mbB}$ in extremes. If $\nu$ has its mass only on the axes, then $\|X\|_{\mbB}$ and $\|Y\|_{\mbB}$ are asymptotically independent, i.e., if one curve shows an extreme behavior, there is negligible probability of the other curve also showing an extreme behavior. If $\nu$ has mass only on the line $\{t(1, 1), t > 0\}$, then $\|X\|_{\mbB}$ and $\|Y\|_{\mbB}$ show asymptotic full dependence, i.e., extreme curves occur simultaneously in $X$ and $Y$.

We remark that the measure $\nu$ on $[0,\infty)^2 \setminus \bzero$ is homogeneous because for any $A \subset [0,\infty)^2 \setminus \{\bzero\}$ and $t>0$,

\begin{align*}
\nu(tA) &= \mu \lbr (x,y) \in \mbB^2_0: (\|x\|_{\mbB},  \|y\|_{\mbB}) \in tA \rbr\\
&= \mu \lbr t(x^{\prime},y^{\prime}) \in \mbB^2_0: (\|x^{\prime}\|_{\mbB},  \|y^{\prime}\|_{\mbB}) \in A \rbr \\
&= t^{-\ag}\mu \lbr (x^{\prime},y^{\prime}) \in \mbB^2_0: (\|x^{\prime}\|_{\mbB},  \|y^{\prime}\|_{\mbB}) \in A \rbr\\
\end{align*}

The joint angular probability measure  $\Gg$ characterizes how the shapes of scaled $X$ and $Y$ are related in extremes. If the extreme curves are exactly proportional, i.e., $X=\la Y$, $\la \neq 0$, the scaled curves share the same extreme functional elements. This means that $\Gg$ concentrates on the “line" $\{(\phi_j,\phi_j),\ j\in \cJ \subset \mbN \}\subset\mbS^2$, where $\{\phi_j,\ j \ge 1\}$ is a set of orthonormal elements in $\mbS$. If the shapes of two curves do not match in extremes, then $\Gg$ concentrates on $\{(\phi_j,\phi_k)\}\subset\mbS^2$, where $ j\in \cJ \subset \mbN$,  $k \in \cK\subset \mbN$, $\cJ \cap \cK = \emptyset$. This situation corresponds to vanishing extremal covariance defined in Section <ref>.

The marginal extreme behavior of $X$ can be obtained by integrating all possible values of $Y$ in def-RVXY, or $\|Y\|_{\mbB}, Y/\|Y\|_{\mbB}$ in pol-RVXY.  Then $X$ has its marginal measure $\mu_X$, and equivalently $X \in RV(-\ag, \Gg_X)$, where $\Gg_X$ is the marginal angular measure of $X$. Similarly, $Y$ has its marginal $\mu_Y$, and equivalently $Y \in RV(-\ag, \Gg_Y)$, where $\Gg_Y$ is the marginal angular measure of $Y$.

§ EXTREMAL CORRELATION COEFFICIENT FOR FUNCTIONAL DATA

The range of $\sg_{XY}$, extremal covariance between $X$ and $Y$, depending on the extremal dependence between $\|X\|$ and $\|Y\|$ and the level of similarity between $X/\|X\|$ and $Y/\|Y\|$ in extremes.

$\|X\|$ and $\|Y\|$     2c$\|X\|$ and $\|Y\|$

Asymptotic Independence     2cAsymptotic Dependence

$X/\|X\|$ and $Y/\|Y\|$

$\sg_{XY}= 0$      look similar     $\sg_{XY} >0$

look orthogonal     $\sg_{XY}\approx 0$

look opposite     $\sg_{XY}<0$

In this section, we introduce the extremal correlation coefficient for functional data. It focuses on the extreme part of the joint distribution of regularly varying random elements $X$ and $Y$ in $L^2$. It measures the tendency of paired curves to exhibit similar extreme patterns by computing a suitable inner product between $X$ and $Y$ conditional on large $[X, Y]^{\top}$.

Given a regularly varying bivariate random element $[X,Y]^{\top}$ in $L^2\times L^2$ with joint exponent measure $\mu$, we define the extremal covariance between $X$ and $Y$ by

\begin{equation} \label{e:cov}
\sg_{XY} = \int_{\|x\| \vee \|y\|>1} \lip x, y \rip  \mu(dx, dy).
\end{equation}

Recall that by mu1, $\mu I_{\{\|x\| \vee \|y\|>1\}}$ is a probability measure. The extremal covariance is thus an extreme analog of the classic covariance in that $\sg_{XY}$ measures how much two random curves vary together in extremes. To see this more closely, we recall the transformation $T$ defined in T and the relation $\mu\circ T^{-1} = \nu \times \Gg$ in nuGg. It then follows from the change of variables that

\begin{equation} \label{e:cov_pol}
\sg_{XY} = \int_{r_X \vee r_Y >1} r_Xr_Y \ \nu (dr_X,dr_Y) \int_{\mbS^2} \lip \thg_X, \thg_Y \rip \Gg (d\thg_X,d\thg_Y).
\end{equation}

The extremal covariance of $X$ and $Y$ can be thus factorized into the extremal dependence between $\|X\|$ and $\|Y\|$ and the level of similarity of the shapes between $X/\|X\|$ and $Y/\|Y\|$. If $\|X\|$ and $\|Y\|$ are asymptotically independent, i.e., extreme curves do not occur simultaneously, then $\nu$ concentrates on the coordinate axes. This implies that $\int_{r_X \vee r_Y >1} r_Xr_Y \ \nu (dr_X,dr_Y)=0$, so $\sg_{XY}=0$ regardless of how $X$ and $Y$ look like in extremes. If $\|X\|$ and $\|Y\|$ are asymptotically dependent, i.e., extreme norms tend to occur simultaneously, then $\int_{r_X \vee r_Y >1} r_Xr_Y \ \nu (dr_X,dr_Y)> 0$ and so there are three possible ranges for $\sg_{XY}$ depending on the relative shape of extreme $X/\|X\|$ and $Y/\|Y\|$: 1) $\sg_{XY}>0$ if the shapes look similar, 2) $\sg_{XY} \approx 0$ if the shapes do not match, i.e., the curves are orthogonal, or 3) $\sg_{XY}<0$ if they look opposite. These properties are summarized in Table <ref>.

We further define the extremal correlation coefficient by

\begin{equation} \label{e:cor}
\rho_{XY} = \frac{\sg_{XY}}{\sg_{X}\sg_{Y}},
\end{equation}

\begin{align*}
\sg_{X} = \lbr \int_{\|x\| \vee \|y\|>1}  \|x\|^2   \mu(dx,dy) \rbr^{1/2}, \ \ \ \
\sg_{Y} = \lbr \int_{\|x\| \vee \|y\|>1}  \|y\|^2   \mu(dx,dy) \rbr^{1/2}.
\end{align*}

The coefficient $\rho_{XY}$ has  properties  analogous to the classic correlation coefficient: 1) $-1\le \rho_{XY} \le 1$, 2) $\rho_{XY}$ measures the strength and direction of linear relationships between $X$ and $Y$ in extremes, 3) If $X$ and $Y$ are independent, then $\rho_{XY}=0$ since independence implies asymptotic independence between $\|X\|$ and $\|Y\|$.

To motivate our estimation approach, we first show that  $\sg_{XY}$ is a  limit of the expected inner product of  $X$ and $Y$ conditional on large values of $[X,Y]^{\top}$.

Let $[X,Y]^{\top}$ be a regularly varying random element in $L^2 \times L^2$. Then,

\[
\sg_{XY} = \lim_{n\to \infty}E\lb \lip \frac{X}{b(n)}, \frac{Y}{b(n)} \rip \Bigg|\|X\| \vee \|Y\| >b(n)\rb.
\]

Considering $f: L^2 \times L^2 \to \mbR$ defined by $(x, y) \to \lip x, y \rip I_{\|x\| \vee \|y\|>1}$, we have that

\begin{align*}
&E\lb \lip {b(n)}^{-1}X, {b(n)}^{-1}Y \rip | \|X\| \vee \|Y\| >{b(n)}\rb \\
&=\frac{1}{P(\|X\|\vee\|Y\|>{b(n)})} E\lb \lip {b(n)}^{-1}X, {b(n)}^{-1}Y \rip I_{\|X\| \vee \|Y\| >{b(n)}}\rb\\
&=\int_{L^2 \times L^2} f(x,y) \frac{P({b(n)}^{-1}X \in dx, {b(n)}^{-1}Y \in dy)}{P(\|X\|\vee\|Y\|>{b(n)})}.
\end{align*}

Note that $f$ is bounded and vanishes on a neighborhood of $\bzero$ in $L^2 \times L^2$. Also, the discontinuity set of $f$ is the boundary of $\mathcal A_1 = \{(x,y) \in L^2 \times L^2: \|x\| \vee \|y\| \ge 1\}$, and it follows from  Lemma <ref> that $\mu(\partial \mathcal A_1)=0$. Therefore, by def-RVXY and Lemma A.1 of [27], we get the claim.

Based on Proposition <ref>, we propose an estimator for $\sg_{XY}$ defined by

\begin{equation} \label{e:ec}
\hat{\sg}_{n,k} =\frac{1}{k}\sum_{i=1}^n \lip \frac{X_i}{R_{(k)}}, \frac{Y_i}{R_{(k)}} \rip I_{R_i \ge R_{(k)} },
\end{equation}

where $[X_i, Y_i]^{\top}, 1 \le i \le n$, are i.i.d. copies of $[X,Y]^{\top}$,   $R_i:=\|X_i\| \vee \|Y_i\|$ and $R_{(k)}$ is the $k$th largest order statistic with the convention $R_{(1)} = \max \{ R_1, \ldots, R_n \}$. An estimator for $\rho_{XY}$ is then defined by

\begin{equation} \label{e:ecc}
\hat{\rho}_{n,k} = \frac{\sum_{i=1}^n \lip X_i, Y_i \rip}{\sqrt{\sum_{i=1}^n \|X_i\|^2} \sqrt{\sum_{i=1}^n \|Y_i\|^2}} \ I_{R_i \ge R_{(k)} }.
\end{equation}

These estimators take only the $k$ largest pairs of $[X_i, Y_i]^{\top}$, $1\le i\le n$, as inputs. This approach falls into so-called peaks-over-threshold framework in that it relies only on $k$ largest observations whose magnitude exceeds a certain threshold. Asymptotic properties in this framework are typically derived as $k$ goes to infinity with $n$, in such a way that $k/n \to 0$. We assume throughout the paper that this condition holds.

We will work under the following assumption.

The bivariate random element $[X,Y]^{\top}$ in  $L^2 \times L^2$ has mean zero and is regularly varying with index $-\ag$, $\ag >2$. The observations $[X_1,Y_1]^{\top}, [X_2,Y_2]^{\top}, \ldots$ are independent copies of $[X,Y]^{\top}$.

We state in the following theorem that the estimator $\hat{\sg}_{n,k}$ is consistent for the extremal covariance. All proofs of the theoretical results introduced in this section are presented in Sections <ref> and <ref> of the Appendix, as they require a number of preliminary  results and technical arguments.

Under Assumption <ref>,

\[
\hat{\sg}_{n,k} \convP \sg_{XY},
\]

where $\hat{\sg}_{n,k}$ and $\sg_{XY}$ are defined in ec and cov, respectively.

From Theorem <ref>, the consistency of $\hat{\rho}_{n,k}$ for $\rho_{XY}$ follows from  Slutsky's theorem.

Under Assumption <ref>,

\[
\hat{\rho}_{n,k} \convP \rho_{XY},
\]

where $\hat{\rho}_{n,k}$ and $\rho_{XY}$  are defined in ecc and cor, respectively.

We end this section with a discussion on the condition $\ag>2$ in Assumption <ref>. The condition is needed because the definition of extremal correlation coefficient presumes the existence of the second moment of the underlying processes. It  can be lifted for the following alternative measure

\begin{equation} \label{e:ga}
\ga_{XY} := \int_{\mbS^2} \lip \thg_X, \thg_Y \rip \Gg(d\thg_X, d\thg_Y),
\end{equation}

which is the angular density factor  in cov_pol. This measure itself could be considered a measure for extremal correlation since $-1 \le \ga_{XY} \le 1$.
Just as Proposition <ref>, it can be shown that

\[
\ga_{XY} = \lim_{n\to \infty} E \lb \lip \frac{X}{\|X\|}, \frac{Y}{\|Y\|} \rip \Bigg| \|X\| \vee \|Y\| >b(n)\rb.
\]

From this relation, an estimator for $\ga_{XY}$ can be defined by

\[
\hat{\ga}_{n,k} = \frac{1}{k} \sum_{i=1}^N  \lip \frac{X_i}{\|X_i\|}, \frac{Y_i}{\|Y_i\|} \rip I_{R_i \ge R_{(k)}},
\]

and its consistency can be proven in almost the same manner as  the proof of Corollary 4.2 of [2],  where $\ag>0$ is assumed.

Although the assumption $\ag>2$ could be relaxed by employing $\ga_{XY}$ as an extremal correlation measure, it should be noted that $\ga_{XY}$ does not account for whether extreme curves $X$ and $Y$ occur simultaneously,  as $\rho_{XY}$ in cor does. To see this, suppose that $X=Z_1\phi$ and $Y=Z_2\phi$, where $Z_1$, $Z_2$ are independent random variables satisfying $P(Z_{1}>z) = z^{-\ag}$, $z>0$, and $\phi$ is any unit norm
element of $L^2$. Then,  $\ga_{XY} = 1$, but $\rho_{XY}=0$.
One can argue that the measure $\rho_{XY}=0$ is more appropriate
because the extremes of $X$ and $Y$ occur independently,
and the objective of our coefficient is to measure similarity
of shapes in paired extreme curves. The identity $\rho_{XY}=0$
holds because the measure  $\nu$ in cov_pol concentrates on the coordinate axes, so

\[
\sg_{XY} = \int_{r_X \vee r_Y >1} r_Xr_Y \ \nu (dr_X,dr_Y) \int_{\mbS^2} \lip \thg_X, \thg_Y \rip \Gg (d\thg_X,d\thg_Y) = 0\times 1=0,
\]

and the condition $\ag >2$ is needed for the first integral to exist.

§ A SIMULATION STUDY

We perform a simulation study to demonstrate that the proposed estimator, $\hat{\rho}_{n,k}$, consistently estimates the  extremal correlation. We generate functional observations in such a way that the theoretical value of $\rho_{XY}$ can be computed analytically, so that we can see how close $\hat{\rho}_{n,k}$ is to the true value.

The design of our study is as follows. Suppose that $Z_1$, $Z_2$ are i.i.d. random variables in $\mbR$ satisfying $P(|Z_{1}|>z) = z^{-\ag}$ with equal chance of $Z_1$ being either negative or positive. Also, let $N_1$, $N_2$, $N_3$ be i.i.d. normal random variables in $\mbR$ with mean 0 and variance 1. Consider $\{\phi_j, j \ge 1 \}$ defined by basis and recall that it is an  orthonormal basis in $L^2([0,1])$. These functions are simulated on a grid of 100 equally–spaced points on the unit interval $[0,1]$. We consider the following data generating processes, for $-1 \le \rho \le 1$,

\begin{align} \label{e:XY}
X(t) &= Z_{1}\phi_1(t) + N_1\phi_2(t) +  N_2\phi_3(t); \\ \nonumber
Y(t) &= \rho Z_{1}\phi_1(t) + (1-\rho^2)^{1/2} Z_{2}\phi_2(t) + N_3\phi_3(t).
\end{align}

This generates extreme curves dominated by the shape of the functional axis $\phi_1$ for $X$ and  by either $\phi_1$ or $\phi_2$ for $Y$.  The following lemma gives an analytic formula for $\rho_{XY}$. Its proof is provided in Section <ref> of the Appendix.

Let $\bZ = [Z_1, Z_2]^{\top}$ be a random vector in $\mbR^2$ consisting of iid components $Z_j$ whose magnitude is regularly varying with tail index $\ag>2$, i.e., for some $c_{+}$, $c_{-}\ge 0$,

\[
P(Z_1 >z) \sim c_{+}z^{-\ag}, \ \ \ \ P(Z_1 <-z) \sim c_{-}z^{-\ag},
\]

where $f(z) \sim g(z)$ if and only if $\lim_{z\to \infty}f(z)/g(z)=1$.
Also, let $\{\phi_j,\ j \ge 1\}$ be a set of orthonormal elements in $\mbS$. Then, for  $X$ and $Y$ in XY,

\[
\rho_{XY} = \frac{\rho}{\{\rho^2 + (1-\rho^2)^{\ag/2}\}^{1/2}}.
\]

Empirical biases (standard errors) of $\hat{\rho}_{n,k}$ in ecc when $\ag=2.1$. The number of upper order statistics is $k = 2\lfloor n^{1/2} \rfloor$.

$\rho_{XY}$     $N=100$      $N=500$     $\rho_{XY}$     $N=100$      $N=500$

0         0.000 (0.051)        0.000 (0.023)

0.1       0.005 (0.086)        0.015 (0.096)     -0.1       -0.008 (0.097)       -0.009 (0.070)

0.2      0.012 (0.138)       0.014 (0.107)     -0.2       -0.009 (0.133)       -0.016 (0.116)

0.3      -0.012 (0.159)       0.012 (0.137)     -0.3        0.002 (0.166)       -0.013 (0.142)

0.4      -0.021 (0.176)       0.007 (0.149)     -0.4        0.010 (0.180)       -0.011 (0.156)

0.5      -0.041 (0.178)       0.002 (0.164)     -0.5        0.030 (0.179)        0.008 (0.161)

0.6      -0.066 (0.180)       -0.011 (0.168)     -0.6        0.056 (0.180)        0.012 (0.162)

0.7      -0.071 (0.179)       -0.029 (0.153)     -0.7        0.076 (0.174)        0.032 (0.160)

0.8      -0.093 (0.155)       -0.040 (0.142)     -0.8        0.097 (0.158)        0.032 (0.138)

0.9      -0.112 (0.127)       -0.043 (0.108)     -0.9        0.111 (0.127)        0.041 (0.102)

1.0      -0.116 (0.066)       -0.040 (0.021)     -1.0        0.118 (0.069)        0.041 (0.021)

We consider $\rho_{XY} \in \{0, \pm.1, \pm.2, \ldots, \pm .9, \pm 1 \}$ and $\ag \in\{2.1, 3,4,5\}$, from which values of $\rho$ can be obtained by Lemma <ref>. For each $\rho$, we generate $[X_i, Y_i]^{\top}$, $1 \le i \le N$, that are i.i.d. copies of $[X,Y]^{\top}$, with sample sizes $N \in \{100, 500\}$. In each case, 1000 replications are generated.

Empirical biases (standard errors) of $\hat{\rho}_{n,k}$ when $\ag=2.1$. The KS method is used to choose optimal $k$s. On average, it selects $k=8$ ($N=100$) and $k=25 \sim 32$ ($N=500$).

$\rho_{XY}$     $N=100$      $N=500$     $\rho_{XY}$     $N=100$      $N=500$

0        -0.001 (0.083)        0.000 (0.041)

0.1       0.026 (0.139)        0.032 (0.127)     -0.1       -0.033 (0.153)       -0.026 (0.117)

0.2       0.042 (0.204)        0.033 (0.156)     -0.2       -0.043 (0.198)       -0.040 (0.179)

0.3       0.022 (0.225)        0.035 (0.190)     -0.3       -0.033 (0.228)       -0.032 (0.195)

0.4       0.013 (0.229)        0.023 (0.193)     -0.4       -0.023 (0.234)       -0.028 (0.204)

0.5      -0.006 (0.232)        0.019 (0.204)     -0.5       -0.006 (0.231)      -0.002 (0.215)

0.6      -0.032 (0.235)        0.006 (0.211)     -0.6        0.021 (0.233)       -0.002 (0.207)

0.7      -0.034 (0.221)       -0.019 (0.188)     -0.7        0.039 (0.217)        0.025 (0.198)

0.8      -0.049 (0.195)       -0.033 (0.173)     -0.8        0.054 (0.192)        0.024 (0.170)

0.9      -0.067 (0.153)       -0.032 (0.122)     -0.9        0.066 (0.146)        0.034 (0.124)

1.0      -0.064 (0.056)       -0.029 (0.026)     -1.0        0.065 (0.059)        0.030 (0.027)

In order to compute $\hat{\rho}_{n,k}$ in ecc, we must select $k$ largest pairs of curves. We first consider $k=2\lfloor N^{1/2} \rfloor$, where $\lfloor x \rfloor$ is the integer part of $x$, to demonstrate the performance of the estimator with a deterministic form of $k$. Additionally, we provide a method for determining the optimal value of $k$ as a guiding tool in practical applications.

Our approach to identifying the optimal value of $k$ is motivated by the theoretical property of the size of pairs, $R_i=\|X_i\| \vee \|Y_i\|$. By Lemma <ref> (i), if $(X, Y)$ are regularly varying in $L^2 \times L^2$ with tail index $\ag$, then $R = \|X\| \vee \|Y\|$ is also regularly varying in $\mbR_+$ with the same tail index $\ag$. Therefore, we choose $k$ that results in successful tail estimation for $R$ in finite samples. In the literature on tail estimation, various methods for selecting upper-order statistics have been introduced, e.g., [11], [10], [8], and [4], just to name a few. We adopt one of the methods proposed by [3]. It chooses $k$ that minimizes the“Kolmogorov-Smirnov"(KS) distance between the empirical tail and the theoretical tail of a Pareto distribution, which is shown in their paper to exhibit relatively good/stable performance in finite samples. This method is implemented by the function mindist of the R package tea.

We now report empirical biases (average minus theoretical value) and standard errors computed as sample standard deviations. For $\ag=2.1$, the results with $k=2\lfloor N^{1/2} \rfloor$ are presented in Table <ref>, and the results with the optimal $k$s based on the KS method in Table <ref>. We put the results with the KS method when $\ag \in \{3,4,5\}$ in Section <ref> of Appendix since they show similar results as when $\ag=2.1$. The conclusions are summarized as follows.

* The estimator is consistent as both the bias and the standard errors decrease with increasing sample sizes, across almost all values of $\rho_{XY}$. The data-driven method for selecting $k$ yields slightly higher standard errors than the deterministic method, but it produces relatively smaller biases when $|\rho_{XY}|$ is close to 1.

* The bias tends to increase in magnitude as $|\rho_{XY}|$ approaches 1. This could be due to the effect of the boundary, $\rho_{XY} \in \{-1, 1\}$. These barriers make the estimator underestimate the true value.

* The standard errors are observed  to be non-uniform across $\rho_{XY}$, they roughly behave like a quadratic function of $\rho_{XY}$ with its peak at $\pm.5$.

The last finding suggests  that the asymptotic variance of $\hat{\rho}_{n,k}$ could be proportional to $|\rho_{XY}|(1-|\rho_{XY}|)$, just like for  the classic correlation coefficient. The derivation of the asymptotic distribution of
$\hat \rho_{XY}$ is postponed to future work.

§ APPLICATIONS TO FINANCIAL AND CLIMATE DATA

In this section, we compute the extremal correlation coefficient
for a number of paired functional data samples that fall into
two categories: intraday returns and daily temperatures.
Our objective is to show that the coefficient provides meaningful and useful information.

§.§ Extremal dependence of  intraday returns
on sector ETFs

In this section, we study pairwise extremal dependence of cumulative intraday return curves (CIDRs) of  Exchange Traded Funds (ETFs) reflecting
performance of key sectors of the U.S. economy. We work with  nine
Standard & Poor's  Depositary Receipt ETFs listed in Table <ref>.
Our objective is to
measure the tendency of paired CIDRs to exhibit similar extreme daily trajectories during the market decline  caused by the COVID-19 pandemic.
The CIDRs are defined as follows. Denote by $P_i(t)$ the price of an asset
on trading day $i$ at time $t$. For the assets in our example, $t$ is time in minutes between 9:30 and 16:00 EST (NYSE opening times) rescaled to the unit interval $(0,1)$. The CIDR on day $i$ is  the curve

\[
R_i(t) = \ln P_i(t) - \ln P_i(0), \ \ \ t\in [0,1],
\]

where $P_i(0)$ is the opening price on day $i$. The curves $R_i$ show how the return accumulates over the trading day, see  Figure <ref>.
We consider all full trading days between  Jan 02, 2020 and  July 31, 2020 ($N=147$).

The nine sector ETFs  and the estimates $\hat{\ag}$ computed by applying the Hill estimator to $\| R_i\|$  with the function mindist of the R package tea.

Ticker     Sector     $\hat{\ag}$

XLY    Consumer Discretionary    3.8

XLP    Consumer Staples    2.6

XLE    Energy    4.2

XLF    Financials    4.0

XLV    Health Care    3.9

XLI    Industrials    3.7

XLB    Materials    3.4

XLK    Technology    4.7

XLU    Utilities    3.8

The CIDR of four ETFs on the four most extreme days during the Covid-19 market decline. Curves of matching color and type
represent curves on the same day; XLF is paired with XLK and
XLY with XLB.

Estimates of the pairwise extremal correlation coefficients of CIDRs across the nine sectors.

Recall that the mathematical framework def-RVXY from
which $\rho_{XY}$ is derived  assumes that the marginal distributions of $X$ and $Y$ are tail equivalent.
Using the estimates $\hat\ag$ in Table <ref>
and a power transformation, we standardize all tails to $\ag=2.5$.
For completeness, we recall this  method. Given $X \in RV(-\ag_X, \Gg_X)$ and $Y \in RV(-\ag_Y, \Gg_Y)$, consider the transformation

\[
g_X(x) = \frac{x}{\|x\|^{1-{\ag_X}/\ag}}, \ \ \ \ g_Y(y) = \frac{y}{\|y\|^{1-{\ag_Y}/\ag}},\ \ \ \ x,y \in L^2,
\]

where $\ag$ is a desired tail index. Applying $g_X$ and $g_Y$ to $X$ and $Y$, respectively, makes them tail equivalent because $P(\|g_X(X)\|>\cdot)$ and $P(\|g_Y(Y)\|>\cdot)$  are regularly varying with $-\ag$. Since this method adjusts only the scale of curves, the transformed curves still retain their original shapes.
After applying  the above  transform, we  select an optimal $k$ for each pair using the KS method described in Section <ref> to compute $\hat{\rho}_{n,k}$.

Figure <ref> shows estimates of the pairwise extremal correlation coefficient across the nine ETF sectors. All pairs exhibit positive extremal correlations ($\hat{\rho}_{n,k} = 0.39 \sim 0.96$), and 56% of the pairs have strong extremal correlations above 0.7. We see  that the CIDRs overall exhibit matching  patterns of cumulative
intraday returns on extreme market volatility days  during the Covid-19 market turbulence. To the first approximation, on such days,  almost all sectors drop together or increase together. However, our coefficient reveals more subtle
information as well. For example,  extreme return curves of
XLF (Financials) are exceptionally strongly correlated with extreme curves for
XLV, XLB and XLK (Health Care,  Materials, Technology), but relatively weakly
correlated with XLU (Utilities). We use this example for illustration and
do not aim at an analysis of the stock market or the economy, but we
note that some findings are interesting. One might expect that the
financial sector (mostly banks) will be strongly affected by the
technology sector (mostly large IT companies like Google or Microsoft)
because such mega corporations  dominate the U.S. stock market.
The similarity of extreme  return curves for XLF and XLK is illustrated in the two left panels of Figure <ref>. One could
also expect that the stocks of banks will be less affected by the
performance of utility companies whose revenues are to a large extent
fixed. But it is less obvious that banks are strongly correlated
with Health Care and Materials sectors.  As another comparison, consider
XLY (Consumer Discretionary) and XLB (Materials) that  show a weak extremal correlation, $\hat{\rho}_{n,k}=0.39$. Their extreme curves exhibit dissimilar patterns, see the right two panels in Figure <ref>.

§.§ Extremal correlation between daily temperature curves

The three locations in the United States: Fort Collins, CO; Colorado Springs, CO; Austin, TX. The pairwise extremal correlation of daily temperature curves between the three locations is evaluated.

In this section, we evaluate the tendency of paired daily temperature curves to exhibit similar extreme patterns across three locations in the United States. The three locations are marked  in Figure <ref>. We focus on the pairwise extremal dependence of those curves during the 2021 heat wave. Although this example  focuses on temperature curves, our tool can be used
for analyzing other curves during extreme weather events;   for example,
daily precipitation patterns or river flows during floods. A correlation of extreme data   during   past events may help with planning a resilient
infrastructure that can better withstand the next extreme weather event.

We use hourly temperature measurements provided by the European Centre for Medium-Range Weather Forecasts (ECMWF). The data are part of their ERA5 (Fifth Generation of ECMWF atmospheric reanalyses) dataset, and represent the temperatures of air at 2 meters above the surface of land, sea or inland waters. We refer to [12] for more details on the ERA5 data.  We partition the hourly data into daily curves, with each day's curve beginning at UTC$+0$, to produce concurrent daily temperature curves across locations in different time zones. We denote the temperature (in Celsius) on day $i$ at hour $t$ and at location $s \in \{ {\rm Fort \ Collins},\ {\rm Colorado \ Springs}, \ {\rm Austin}\}$ by
$X_i(s, t)$, $i=1, \ldots, N$.
Figure <ref> depicts examples of daily temperature curves at the three locations. The data are taken from May 17, 2021 to Aug 31, 2021 ($N$ = 107).

Extreme daily temperature curves (in Celsius) during the 2021 heat wave (local time on the x-axis).  Curves of matching color represent the same days when both Fort Collins and Colorado Springs   experienced extreme patterns simultaneously.

Prior to computing $\hat{\rho}_{n,k}$ for each pair of the three locations, daily curves are centered by the mean function, $\bar{X}_N (s,t) = \frac{1}{N}\sum_{i=1}^N X_i(s, t)$, for each location $s$.
We then compute the tail index estimate $\hat{\ag}$  of the
norms $\| X_i(s,t) -\bar{X}_N (s,t)\|$, for each location $s$, shown in Table <ref>. The results suggest that the marginal distributions of those curves across the three locations are not tail-equivalent. We apply the power transformation method, described in Section <ref>, to get the tail index $\ag=2.5$ at all locations. We then apply the KS method, described in Section <ref>, to the centered curves to select an optimal $k$ for each pair.

Tail index estimates $\hat{\ag}$ and estimates of the pairwise extremal correlation coefficients, $\hat{\rho}_{n,k}$, of daily temperature curves across Fort Collins, CO, Colorado Springs, CO, and Austin, TX.

Location     $\hat{\ag}$        Fort     Colorado      Austin

Collins     Springs

Fort Collins     5.5        1    0.94     0.79

Colorado Springs     3.6       0.94     1    0.74

Austin    4.4       0.79     0.74   1

Table <ref> reports estimates of the pairwise extremal correlation coefficient across the three locations. There are positive and strong extremal correlations among all pairs ($\hat{\rho}_{n,k} = 0.74 \sim 0.94$), suggesting a high degree of association between the daily temperature extreme patterns across the three locations, even between different
climatic regions like the Front Range  foothills and the southern edge of the Great Plains.
We see that the proximity in geographical locations corresponds to greater similarity in extreme patterns, showing that  $\hat{\rho}_{n,k}$ is
a meaningful and useful dependence measure.

We Thank Professor Hong Miao of the Department of Finance and Real
Estate at Colorado State University for preprocessing the financial data used in Section <ref>. We thank Professor Joshua French of
the Department of Mathematical and Statistical Sciences at the University of Colorado Denver
for preprocessing the temperature data used in Section <ref>.

Conflicts of interest:  None declared.

P.K. was partially supported by the United States National Science Foundation
grant DMS–2123761.

Data availability

The raw high frequency financial data used in Section <ref> can be acquired from tickdata.com, as well as from many other providers.
The reanalysis temperature data used in Section <ref>
can be downloaded from ecmwf.int. The R code used in this paper can be found at github.com/veritasmih/ecc.

Supplementary material

Supplementary material is available at Journal of the Royal Statistical Society: Series B online

[1]
Billingsley, P.
Convergence of Probability Measures; Second Edition.
Wiley, New York.

[2]
Clémençon, S., Huet, N. and Sabourin, A.
Regular variation in Hilbert spaces and principal component analysis for functional extremes.
Stochastic Processes and their Applications; forthcoming.

[3]
Danielsson, J., Ergun, L. M., de Haan, L. and de Vries, C. G.
Tail index estimation: Quantile driven threshold selection.
Technical Report. Bank of Canada.

[4]
Danielsson, J., de Haan, L., Peng, L. and de Vries, C.G.
Using a bootstrap method to choose the sample fraction in tail index estimation.
Journal of Multivariate Analysis, 76.

[5]
Davis, R. A. and Mikosch, T.
The extremogram: A correlogram for extreme events.
Bernoulli, 15, 977–1009.

[6]
de Haan, L. and Ferreira, A.
Extreme Value Theory.

[7]
Dombry, C. and Ribatet, M.
Functional regular variations, Pareto processes and peaks over threshold.
Statistics and its Interface, 8, 9–17.

[8]
Drees, H. and Kaufmann, E.
Selecting the optimal sample fraction in univariate extreme value estimation.
Stochastic Processes and their Applications, 75, 149–172.

[9]
Dubin, J. A. and Müller, H. G.
Dynamical correlation for multivariate longitudinal data.
Journal of the American Statistical Association, 100, 872–881.

[10]
Hall, P.
Using the bootstrap to estimate mean squared error and select smoothing parameter in nonparametric problems.
Journal of Multivariate Analysis, 32, 177–203.

[11]
Hall, P. and Welsh, A. H.
Adaptive estimates of parameters of regular variation.
The Annals of Statistics, 13, 331–341.

[12]
Hersbach, H., Bell, B., Berrisford, P., Hirahara, S., Horányi, A., Muñoz-Sabater, J., Nicolas, J., Peubey, C., Radu, R., Schepers, D. et al.
The ERA5 global reanalysis.
Quarterly Journal of the Royal Meteorological Society, 146, 1999–2049.

[13]
Hsing, T. and Eubank, R.
Theoretical Foundations of Functional Data Analysis, with an Introduction to Linear Operators.

[14]
Hult, H. and Lindskog, F.
Regular variation for measures on metric spaces.
Publications de l'Institute Mathématique. Nouvelle Série, 80(94), 121–140.

[15]
Janßen, A., Neblung, S. and Stoev, S.
Tail-dependence, exceedance sets, and metric embeddings.
Extremes, 1–39.

[16]
Kim, M. and Kokoszka, P.
Hill estimator of projections of functional data on principal components.
Statistics, 53, 699–720.

[17]
Kim, M. and Kokoszka, P.
Extremal dependence measure for functional data.
Journal of Multivariate Analysis, 189, 104887.

[18]
Kokoszka, P. and Kulik, R.
Principal component analysis of infinite variance functional data.
Journal of Multivariate Analysis, 193, 105123.

[19]
Kokoszka, P. and Reimherr, M.
Introduction to Functional Data Analysis.
CRC Press.

[20]
Kokoszka, P., Stoev, S. and Xiong, Q.
Principal components analysis of regularly varying functions.
Bernoulli, 25, 3864–3882.

[21]
Kokoszka, P. and Xiong, Q.
Extremes of projections of functional time series on data–driven basis systems.
Extremes, 21, 177–204.

[22]
Larsson, M. and Resnick, S. I.
Extremal dependence measure and extremogram: the regularly varying case.
Extremes, 15, 231–256.

[23]
Ledford, A. W. and Tawn, J. A.
Statistics for near independence in multivariate extreme values.
Biometrika, 83, 169–187.

[24]
Ledford, A. W. and Tawn, J. A.
Modelling dependence within joint tail regions.
Journal of the Royal Statistical Society: Series B (Statistical Methodology), 59, 475–499.

[25]
Ledford, A. W. and Tawn, J. A.
Diagnostics for dependence within time series extremes.
Journal of the Royal Statistical Society: Series B (Statistical Methodology), 65, 521–543.

[26]
McDiarmid, C.
In Probabilistic Methods for Algorithmic Discrete Mathematics, pp. 195–248.

[27]
Meiguet, T.
Heavy Tailed Functional Time Series.
Ph.D. Thesis. Universite catholique de Louvain.

[28]
Resnick, S. I.
Extreme Values, Regular Variation, and Point Processes.

[29]
Resnick, S. I.
The extremal dependence measure and asymptotic independence.
Stochastic models, 20, 205–227.

[30]
Resnick, S. I.
Heavy–Tail Phenomena.

[31]
Yeh, C. K., Rice, G. and Dubin, J. A.
Functional spherical autocorrelation: A robust estimate of the autocorrelation of a functional time series.
Electronic Journal of Statistics, 17, 650–687.

Supplementary Material

§ PRELIMINARY RESULTS

In this section, we put together preliminary results needed to prove Theorem <ref>. These results allow us to streamline the exposition of proofs of the main result. Recall that, c.f.,  A, $\mathcal A_r =\{(x,y)\in \mbB_0^2: \|(x,y)\|_{\mbB^2}  \ge r\}$, $r>0$, where $\|(x,y)\|_{\mbB^2} = \|x\|_{\mbB} \vee \|y\|_{\mbB}$.

Suppose $\mu$ is a measure in $M_0(\mbB^2)$ satisfying $\mu(t\cdot)=t^{-\ag}\mu(\cdot)$, $t>0$.  Then, $\mathcal A_r$ is a $\mu$–continuity set, i.e., $\mu(\partial \mathcal A_r)=0$.

We assume $\mu(\partial{\mathcal A_r})>0$ and get a contradiction. Since ${\mathcal A_r}\supset \bigcup_{n\ge 1} \partial(n^{1/\ag}{\mathcal A_r})$, it follows from the homogeneity property of $\mu$ that

\[
\mu({\mathcal A_r}) \ge \sum_{n=1}^{\infty} \mu(\partial(n^{1/\ag}{\mathcal A_r})) = \sum_{n=1}^{\infty} \mu(n^{1/\ag}\partial{\mathcal A_r}) =\sum_{n=1}^{\infty} n^{-1}\mu(\partial{\mathcal A_r})=\infty.
\]

It contradicts to the fact that $\mu$ is boundedly finite.

Recall that $R= \|(X, Y)\|$, $R_i = \|(X_i, Y_i)\|$, and $R_{(k)}$ is the $k$th largest order statistic with the convention $R_{(1)} = \max \{ R_1, \ldots, R_n \}$. Let $b(n)$ be the quantile function such that $P(R>b(n))=n^{-1}$. We claim in the following lemma that regular variation of $[X,Y]^{\top}$ in $L^2\times L^2$ implies regular variation of $R$ in $[0,\infty)$.

Let $M_{+}(0,\infty]$ be the space of Radon measures on $(0,\infty]$, and $\nu_{\ag}(r, \infty] = r^{-\ag}$. If $[X,Y]^{\top}$ is regularly varying in $L^2\times L^2$ according to Definition <ref>, then

(i) $R$ is a nonnegative random variable whose distribution has a regularly varying tail with index $-\ag$,

(ii) $\frac{1}{k} \sum_{i=1}^n I_{R_i / b(n/k)} \convP \nu_{\ag}$, in $M_{+}(0,\infty]$,

(iii) $R_{(k)}/b(n/k) \convP 1$, in $[0,\infty)$,

(iv) $\frac{1}{k} \sum_{i=1}^n I_{R_i / R_{(k)}} \convP \nu_{\ag}$ in $M_{+}(0,\infty]$.

For statement $(i)$, observe that for $r>0$,

\[
\frac{n}{k}P \lp R \ge rb(n/k) \rp = \frac{n}{k}P \lp \frac{(X,Y)}{b(n/k)} \in \cA_{r} \rp,
\]

where $\cA_r$ is defined in A. It then follows from def-RVXY and mu1 that

\[
\frac{n}{k}P \lp R \ge rb(n/k) \rp \to \mu(\cA_r)= r^{-\ag}\mu(\cA_1)=r^{-\ag}.
\]

Therefore, by Theorem 3.6 of [30], $(i)$ holds. By Theorem 4.1 of [30], it can be shown that $(i)$ implies $(ii)$. Also, it follows from Step 1 in the proof of Theorem 4.2 of [30] that $(ii)$ implies $(iii)$ and from Step 2 in that proof that  $(ii)$ and $(iii)$ imply $(iv)$.

The following lemma is used to prove Lemmas <ref> and <ref>.

Suppose $\ga_n$ converges vaguely to $\nu_\ag$ in $M_+(0,\infty]$. Then for any compact interval $K\subset (0,\infty]$,

\[
\int_K r^2\ga_n(dr) \to \int_K r^2 \nu_\ag(dr).
\]

Since the function $r\mapsto r^2 I_{K}$ is not continuous,
we use an approximation argument. Set $K= [a, b]$, for $0<a<b \le \infty$. Construct
compact intervals $K_j \searrow K$ and nonnegative continuous functions
$f_j$ such that $I_K \le f_j \le I_{K_j}$. By the triangle inequality,

\begin{align*}
\left | \int_K r^2 \ga_n(dr) - \int_K r^2 \nu_\ag(dr)\right |
&\le  \left | \int r^2 I_K(r) \ga_n(dr) - \int r^2 f_j(r) \ga_n (dr)\right |\\
& \ \ +
\left | \int r^2 f_j(r) \ga_n (dr)  - \int r^2 f_j(r) \nu_\ag  (dr)\right |\\
& \ \ +
\left | \int r^2 f_j(r) \nu_\ag  (dr)  - \int r^2 I_K(r) \nu_\ag(dr)\right |\\
&=: A_{n, j}^{(1)} + A_{n, j}^{(2)} + A_{j}^{(3)}.
\end{align*}

Fix $\ep > 0$. There is $j^\star$ such that for $j \ge j^\star$,

\[
A_{j}^{(3)} \le c\int \lb f_j(r) - I_K(r) \rb \nu_\ag(dr)
\le c \nu_\ag( K_j \setminus K^\circ) < \ep/2,
\]

where $c=b^2I_{b\neq \infty}+a^2I_{b=\infty}$.
Similarly $ A_{n, j}^{(1)} \le c \ga_n ( K_j \setminus K^\circ)$,
so for every fixed $j$,

\[
\limsup_{n\to \infty} A_{n, j}^{(1)} \le M^2
\limsup_{n\to \infty} \ga_n( K_j \setminus K^\circ)
\le M^2 \nu_\ag( K_j \setminus K^\circ)
\]

because $K_j \setminus K^\circ$ is compact, cf. Proposition 3.12 in
[28]. Thus,

\[
\limsup_{n\to\infty}
\left | \int_K r^2 \ga_n(dr) - \int_K r^2 \nu_\ag(dr)\right |
\le  \ep + \limsup_{n\to\infty} A_{n, j^\star}^{(2)} = \ep.
\]

Since $\ep$ is arbitrary, we get the claim.

The following two lemmas are used to prove Lemma <ref> and Proposition <ref>.

Under Assumption <ref>, for any $M>0$,

\[
\frac{n}{k} E \lb \lp \frac{R}{b(n/k)} \rp^2 I_{R \ge Mb(n/k) }  \rb  \to \frac{\ag}{\ag-2}M^{2-\ag}.
\]

Observe that

\[
\frac{n}{k} E \lb \lp \frac{R}{b(n/k)} \rp^2 I_{R \ge Mb(n/k) }  \rb
=\int_M^{\infty} r^2  \frac{n}{k}P\lp \frac{R}{b(n/k)} \in dr\rp,
\]

\[
\frac{\ag}{\ag-2}M^{2-\ag} = \int_{M}^{\infty} r^2 \nu_{\ag}(dr).
\]

By Lemma <ref> (i), we have that in $M_+(0,\infty]$

\[
\frac{n}{k}P\lp \frac{R}{b(n/k)} \in \cdot \rp \convv \nu_{\ag}.
\]

Therefore, we get the claim by Lemma <ref> with $K= [M, \infty]$.

The function $h$ on $M_+(0,\infty]$ defined by $h(\ga) = \int_1^M r^2 \ga(dr)$ is continuous at $\nu_\ag$.

Suppose $\ga_n$ converges vaguely to $\nu_\ag$. Then, by Lemma <ref> with $K= [1, M]$, it can be shown that

\[
\lim_{n\to\infty} \int_1^M r^2 \ga_n(dr) = \int_1^M r^2 \nu_\ag(dr).
\]

The following lemma is the key argument to prove Proposition <ref>.

Under Assumption <ref>, the following statements hold:

\begin{align}
\frac{1}{k}\sum_{i=1}^n \lp \frac{R_i}{R_{(k)}} \rp^{2}   I_{R_i \ge R_{(k)} }  &\convP \frac{\ag}{\ag-2}; \label{e:2k}\\
\frac{1}{k}\sum_{i=1}^n \lp \frac{R_i}{b(n/k)} \rp^{2}   I_{R_i \ge b(n/k) } &\convP \frac{\ag}{\ag-2}. \label{e:2b}
\end{align}

The proofs for 2k and 2b are almost the same, so we only prove 2k to save space. Let $\hat{\ga}_{n,k}=\frac{1}{k} \sum_{i=1}^n I_{R_i / R_{(k)}}$, and recall that $\hat{\ga}_{n,k} \convP \nu_{\ag}$ (see Lemma <ref> (iv)). Since

\[
\frac{1}{k}\sum_{i=1}^n \lp \frac{R_i}{R_{(k)}} \rp^{2}   I_{R_i \ge R_{(k)} }= \int_1^\infty r^2 \hat{\ga}_{n,k}(dr),
\]

we need to show that

\[
\int_1^\infty r^2  \hat{\ga}_{n,k}(dr)
\convP \int_1^\infty r^{2} \nu_\ag(dr)=\frac{\ag}{\ag-2}.
\]

To prove this convergence, we use the second converging together theorem, Theorem 3.5
in [30], (also stated as Theorem 3.2 of [1]). Let

\begin{align*}
V_{n,k} = \int_1^\infty r^2 \hat{\ga}_{n,k}(dr),  \ \ \
&  V= \int_1^\infty r^2 \nu_\ag(dr);\\
V_{n,k}^{(M)} = \int_1^M r^2 \hat{\ga}_{n,k}(dr),  \ \ \
& V^{(M)}= \int_1^M r^2 \nu_\ag(dr).
\end{align*}

To show the desired convergence $V_{n,k} \convP V$ (equivalently, $V_{n,k} \convd V$), we must verify that

\begin{equation} \label{e:VnM}
\forall \ M > 1, \ \ \ V_{n,k}^{(M)} \convd V^{(M)}, \ \ \ \ {\rm as} \ n \to\infty;
\end{equation}

\begin{equation} \label{e:VM}
V^{(M)} \convd V, \ \ \ {\rm as} \ M\to\infty;
\end{equation}

\begin{equation} \label{e:bill}
\forall \ \eg > 0, \ \ \ \lim_{M\to\infty} \limsup_{n\to\infty}
P\lp | V_{n,k}^{(M)}- V_{n,k} | > \eg \rp = 0.
\end{equation}

Convergence VnM follows from Lemma <ref> (iv)
and Lemma <ref>. Convergence VM holds since for $\ag>2$

\[
\int_M^\infty r^2 \nu_\ag(dr)=\int_M^\infty r^2 \ag r^{-\ag-1} dr =\frac{\ag}{\ag-2} M^{2-\ag} \to 0, \ \ \ \ {\rm as}\ M \to \infty.
\]

It remains to show that $\forall \eg>0$,

\[
\lim_{M\to\infty} \limsup_{n\to\infty}
P\lp | V_{n,k}^{(M)}- V_{n,k} | > \eg \rp =  \lim_{M\to\infty} \limsup_{n\to\infty}
P\lp \int_M^\infty r^2 \hat\ga_{n,k}(dr) > \eg \rp = 0.
\]

Fix $\eg> 0$ and $\eta> 0$. Observe that

\[
P\lp \int_M^\infty r^2 \hat\ga_{n,k}(dr) > \eg \rp
\le Q_1(n) + Q_2(n),
\]

\[
Q_1(n) = P\lp \int_M^\infty r^2 \hat\ga_{n,k}(dr) > \eg, \
\left | \frac{R_{(k)}}{b(n/k)} - 1 \right |  < \eta \rp,\ \
Q_2(n) = P \lp \left | \frac{R_{(k)}}{b(n/k)} - 1 \right |  \ge  \eta\rp.
\]

By Lemma <ref> (iii), $ \limsup_{n\to\infty} Q_2(n) = 0$. For $Q_1(n)$, we start with the bound

\begin{align*}
&\le P\lp  \int_M^\infty r^2 \hat\ga_{n,k}(dr) > \eg, \
\frac{R_{(k)}}{b(n/k)} > 1- \eta \rp\\
&= P \lp  \int_M^\infty r^2 \frac{1}{k} \sum_{i=1}^n
I_{R_i/ R_{(k)} \in dr}  > \eg, \
\frac{R_{(k)}}{b(n/k)} > 1- \eta \rp.
\end{align*}

Conditions $R_i/ R_{(k)} > M$ and ${R_{(k)}}/{b(n/k)} > 1- \eta$
imply ${R_{i}}/{b(n/k)} > M(1- \eta)$, so

\begin{align*}
&\le P \lp  \int_{M(1-\eta)}^\infty r^2 \frac{1}{k} \sum_{i=1}^n
I_{R_i/b(n/k) \in dr}   > \eg \rp\\
&= P \lp  \frac{1}{k} \sum_{i=1}^n  \lp \frac{R_i}{b(n/k)} \rp^2
I_{R_i \ge M(1-\eta)b(n/k) } >\eg\rp.
\end{align*}

Then, it follows from Markov's inequality and Lemma <ref> that

\[
Q_1(n) \le \frac{1}{\eg} \frac{n}{k} E \lb \lp \frac{R_1}{b(n/k)} \rp^2
I_{R_1 \ge M(1-\eta)b(n/k) } \rb  \to \frac{1}{\eg} \frac{\ag}{\ag-2} \{M(1-\eta)\}^{2-\ag}, \ \ \ {\rm as} \ n \to \infty.
\]

This bound goes to 0 as $M \to \infty$ since $\ag>2$.

The following lemma follows from Theorem 3.8 of [26]. It states a Bernstein type inequality, which is the key technique to prove Proposition <ref>.

Let $\bZ_n=(Z_1, \ldots, Z_n)$ with the $Z_i$ taking values in a
Lebesgue measurable subset $\cZ$ of an Euclidean space. Let $f$ be a real-valued function defined on $\cZ^n$. For $(z_1, \ldots, z_i) \in \cZ^{i}$, $1 \le i \le n$, put

\begin{equation} \label{e:gi}
g_i(z_1, \ldots, z_i):=E\lb f(\bZ_n)| Z_j=z_j, 1 \le j \le i\rb -E\lb f(\bZ_n)| Z_j=z_j, 1 \le j \le i-1\rb.
\end{equation}

Define the maximum deviation by

\begin{equation} \label{e:b}
b := \max_{1\le i\le n}\sup_{(z_1, \ldots, z_i) \in \cZ^i} g_i(z_1, \ldots, z_i),
\end{equation}

and define the supremum sum of variances by

\begin{equation} \label{e:v}
\hat{v} := \sup_{(z_1, \ldots, z_n) \in \cZ^n} \sum_{i=1}^n \var\lb g_i(z_1, \ldots, z_{i-1}, Z_i^{\prime})\rb,
\end{equation}

where $Z_i^{\prime}$ is an independent copy of $Z_i$ conditional on $Z_j=z_j$, $1 \le j \le i-1$. If $b$ and $\hat{v}$ are finite, then for any $\eg \ge 0$,

\[
P \lp f(\bZ_n) - E[f(\bZ_n)] \ge t \rp \le \exp \lp \frac{-\eg^2}{2(\hat{v}+b\eg/3)}  \rp.
\]

§ PROOF OF THEOREM <REF> IN SECTION <REF>

Recall ec, i.e., the definition:

\[
\hat{\sg}_{n,k} =\frac{1}{k}\sum_{i=1}^n \lip \frac{X_i}{R_{(k)}}, \frac{Y_i}{R_{(k)}} \rip I_{R_i \ge R_{(k)} }.
\]

To prove the consistency of $\hat{\sg}_{n,k}$ for the extremal covariance $\sg_{XY}$, we consider the following sequence of random variables

\begin{align} \label{e:sg}
\sg_{n,k} &:=\frac{1}{k}\sum_{i=1}^n \lip \frac{X_i}{b(n/k)}, \frac{Y_i}{b(n/k)} \rip I_{R_i \ge b(n/k) }.
\end{align}

Note that $\sg_{n,k}$ is not observable since $b(\cdot)$ is unknown. However, $b(n/k)$ can be estimated by its consistent estimator $R_{(k)}$, and it can be shown that replacing $b(n/k)$ by $R_{(k)}$ ensures that the difference between $\sg_{n,k}$ and $\hat{\sg}_{n,k}$ is asymptotically negligible, which will be shown in Proposition <ref>. Thus, the key argument for establishing the consistency is to show that $\sg_{n,k}$ converges in probability to $\sg_{XY}$, which is proven in the following proposition.

Under Assumption <ref>,

\[
\sg_{n,k} \convP \sg_{XY}.
\]

\begin{equation} \label{e:sgbar}
\bar{\sg}_{n,k} := E\lb \lip \frac{X_1}{b(n/k)}, \frac{Y_1}{b(n/k)} \rip \Bigg|\|X_1\| \vee \|Y_1\| >b(n/k)\rb.
\end{equation}

Then, by Proposition <ref>, $\bar{\sg}_{n,k} \to \sg_{XY}$, so it remains to show that $|\sg_{n,k} - \bar{\sg}_{n,k}| \convP 0$.

Let $\bZ_n = (Z_1, \ldots, Z_n)$, where $Z_i = (X_i, Y_i)$, and $\bz_n = (z_1, \ldots, z_n)$, where $z_i = (x_i, y_i)$, for $1 \le i \le n$. Consider a map $f: (L^2 \times L^2)^n \to \mbR$ defined by

\[
f(\bz_n) := \lmo \frac{1}{k}\sum_{i=1}^n \lip \frac{x_i}{b(n/k)}, \frac{y_i}{b(n/k)} \rip I_{r_i \ge b(n/k) }- \frac{n}{k}E\lb \lip \frac{X_1}{b(n/k)}, \frac{Y_1}{b(n/k)} \rip I_{R_1 >b(n/k)}\rb \rmo.
\]

Then, we have that

\[
|\sg_{n,k} - \bar{\sg}_{n,k}| = f(\bZ_n) - E[f(\bZ_n)] + E[f(\bZ_n)].
\]

We aim to show that $f(\bZ_n) - E[f(\bZ_n)] \convP 0$ and $E[f(\bZ_n)] \to 0$.

To establish the convergence, $f(\bZ_n) - E[f(\bZ_n)] \convP 0$, we use the Bernstein type concentration inequality in Lemma <ref>. Since the $(X_i,Y_i)$ are independent, the deviation function in gi has the following form

\[
g_i(z_1, \ldots, z_i)= E \lb f(z_1, \ldots, z_{i-1},z_i, Z_{i+1}, \ldots, Z_n)-f(z_1, \ldots, z_{i-1},Z_i, Z_{i+1}, \ldots, Z_n) \rb.
\]

Then, using the fact that $\lmo |x|- |y| \rmo \le |x-y|$, we have that

\begin{align*}
g_i(z_1, \ldots, z_i)
&\le \frac{1}{k} E\lb \lmo  \lip \frac{x_i}{b(n/k)}, \frac{y_i}{b(n/k)} \rip I_{r_i \ge b(n/k) } - \lip \frac{X_i}{b(n/k)}, \frac{Y_i}{b(n/k)} \rip I_{R_i \ge b(n/k) }       \rmo \rb \\
&\le \frac{1}{k} \lbr \frac{|\lip x_i, y_i\rip|}{b(n/k)^2} + \frac{k}{n} \frac{n}{k} E \lb \lp \frac{R_i}{b(n/k)} \rp^2 I_{R_i \ge b(n/k) }  \rb \rbr\\
&\le \frac{1}{k} \lbr \frac{\|x_i\| \| y_i\|}{b(n/k)^2} + \frac{n}{k} E \lb \lp \frac{R_i}{b(n/k)} \rp^2 I_{R_i \ge b(n/k) }  \rb \rbr.
\end{align*}

Since $(x_i,y_i) \in L^2\times L^2$ and $\frac{n}{k} E \lb \lp {R_i}/{b(n/k)} \rp^2 I_{R_i \ge b(n/k) }  \rb \to \ag/(\ag-2)$ by Lemma <ref>,  we have that $g_{i}(z_1, \ldots, z_i) \le {c_1}/{k}$, for some constant $c_1>0$. Therefore, the maximum deviation $b$ in b is bounded by $c_1/k$.

Next we investigate the upper bound for the sum of variances $\hat{v}$ in v. Since $E[g_i(z_1, \ldots, z_{i-1}, Z_i^{\prime})]=0$ by the law of total probability, we have that

\begin{align*}
&\var\lb g_i(z_1, \ldots, z_{i-1}, Z_i^{\prime})\rb \\
&= E [g_i^2(z_1, \ldots, z_{i-1}, Z_i^{\prime})]\\
&= E \lb \lbr  f(z_1, \ldots, z_{i-1},Z_i^{\prime}, Z_{i+1}, \ldots, Z_n)-f(z_1, \ldots, z_{i-1},Z_i, Z_{i+1}, \ldots, Z_n) \rbr^2 \rb \\
&\le \frac{1}{k^2} E\lb \lbr  \lip \frac{X_i^{\prime}}{b(n/k)}, \frac{Y_i^{\prime}}{b(n/k)} \rip I_{R_i^{\prime} \ge b(n/k) } - \lip \frac{X_i}{b(n/k)}, \frac{Y_i}{b(n/k)} \rip I_{R_i \ge b(n/k) }       \rbr^2 \rb \\
&\le \frac{2}{k^2} E \lb \lip \frac{X_i}{b(n/k)}, \frac{Y_i}{b(n/k)} \rip^2 I_{R_i \ge b(n/k) }   \rb \\
&\le \frac{2}{k^2} \lbr  \frac{k}{n} \frac{n}{k} E \lb  \lp \frac{R_i}{b(n/k)} \rp^2  I_{R_i \ge b(n/k) }  \rb \rbr.
\end{align*}

It then again follows from Lemma <ref> that $\var\lb g_i(z_1, \ldots, z_{i-1}, Z_i^{\prime})\rb  \le c_2/(nk)$ for some $c_2>0$. Then the supremum sum of variances $\hat{v}$ is bounded above by $c_2/k$. Therefore by Lemma <ref>, for any $\eg>0$

\[
P \lp f(\bZ_n) - E[f(\bZ_n)]  \ge \eg \rp \le \exp \lp \frac{-k\eg^2}{c_1+c_2\eg/3} \rp.
\]

If we apply this inequality to $-f(\bZ_n)$, then we obtain the following `two-sided' inequality

\[
P \lp |f(\bZ_n) - E[f(\bZ_n)]|  \ge \eg \rp \le 2\exp \lp \frac{-k\eg^2}{c_1+c_2\eg/3} \rp.
\]

From this, we obtain that  $f(\bZ_n) - E[f(\bZ_n)] \convP 0$.

Next, to show $E[f(\bZ_n)] \to 0$, we set, for $1 \le i \le n$

\[
\Delta_i =  \lip \frac{X_i}{b(n/k)}, \frac{Y_i}{b(n/k)} \rip I_{R_i \ge b(n/k) }- E\lb \lip \frac{X_1}{b(n/k)}, \frac{Y_1}{b(n/k)} \rip I_{R_1 >b(n/k)}\rb.
\]

Then, we have that

\begin{align*}
E \lb f(\bZ_n) \rb = \frac{n}{k} E \lb \lmo \frac{1}{n} \sum_{i=1}^n \Delta_i \rmo \rb &\le \frac{n}{k} \lbr E \lb \lp \frac{1}{n} \sum_{i=1}^n \Delta_i \rp^2 \rb \rbr^{1/2}\\
&= \frac{n}{k} \lbr E \lb   \frac{1}{n^2} \sum_{i=1}^n \Delta_i^2 + \frac{1}{n^2} \sum_{i \neq j} \Delta_i\Delta_j \rb \rbr^{1/2}.
\end{align*}

Since the $\Delta_i$ are independent, $E[\Delta_i\Delta_j]=0$, for $i \neq j$. Therefore,

\begin{align*}
&E \lb f(\bZ_n) \rb \\
&\le \frac{\sqrt{n}}{k} \lbr E \lb    \Delta_1^2  \rb \rbr^{1/2} \\
&= \frac{\sqrt{n}}{k} \lbr  E  \lb \lp \lip \frac{X_1}{b(n/k)}, \frac{Y_1}{b(n/k)} \rip I_{R_1 \ge b(n/k) } - E\lb \lip \frac{X_1}{b(n/k)}, \frac{Y_1}{b(n/k)} \rip I_{R_1 >b(n/k)}\rb  \rp^2\rb \rbr^{1/2}\\
& = \frac{\sqrt{n}}{k} \lbr  \var \lb  \lip \frac{X_1}{b(n/k)}, \frac{Y_1}{b(n/k)} \rip I_{R_1 \ge b(n/k) } \rb \rbr^{1/2} \\
& \le \frac{\sqrt{n}}{k} \lbr  E  \lb \lip \frac{X_1}{b(n/k)}, \frac{Y_1}{b(n/k)} \rip^2 I_{R_1 \ge b(n/k) } \rb \rbr^{1/2}\\
& \le  \frac{\sqrt{n}}{k}  \lbr E \lb  \lp \frac{R_1}{b(n/k)} \rp^2  I_{R_1 \ge b(n/k) }  \rb \rbr^{1/2}.
\end{align*}

Therefore, by Lemma <ref> we have that

\[
E \lb f(\bZ_n) \rb \le \frac{\sqrt{n}}{k}  \lbr  \frac{k}{n} \frac{n}{k} E \lb  \lp \frac{R_1}{b(n/k)} \rp^2  I_{R_1 \ge b(n/k) }  \rb \rbr^{1/2} \le \frac{c_3}{\sqrt{k}},
\]

for some $c_3>0$, which completes the proof.

Under Assumption <ref>,

\[
|\hat{\sg}_{n,k} - \sg_{n,k} | \convP 0.
\]

Consider the following decomposition

\[
|\hat{\sg}_{n,k} - \sg_{n,k} | \le P_1(n) + P_2(n),
\]

\begin{align*}
&P_1(n) :=\lmo \frac{1}{k}\sum_{i=1}^n \lip \frac{X_i}{R_{(k)}}, \frac{Y_i}{R_{(k)}} \rip  \lbr I_{R_i \ge R_{(k)} } - I_{R_i \ge b(n/k) }\rbr \rmo,\\
&P_2(n) :=\lmo \frac{1}{k}\sum_{i=1}^n \lbr \lip \frac{X_i}{R_{(k)}}, \frac{Y_i}{R_{(k)}} \rip - \lip \frac{X_i}{b(n/k)}, \frac{Y_i}{b(n/k)} \rip \rbr
I_{R_i \ge b(n/k) } \rmo.
\end{align*}

We will show that each of the two parts goes to 0. We first focus on $P_1(n)$. Observe that

\begin{align*}
P_1(n) &\le \lp\frac{b(n/k)}{R_{(k)}} \rp^2  \frac{1}{k}\sum_{i=1}^n  \lmo \lip \frac{X_i}{R_i}, \frac{Y_i}{R_i} \rip \rmo \lp \frac{R_i}{b(n/k)} \rp^{2} \lmo  I_{R_i \ge R_{(k)} } - I_{R_i \ge b(n/k) } \rmo \\
&\le \lp\frac{b(n/k)}{R_{(k)}} \rp^2  \frac{1}{k}\sum_{i=1}^n  \lp \frac{R_i}{b(n/k)} \rp^{2} \lmo  I_{R_i \ge R_{(k)} } - I_{R_i \ge b(n/k) } \rmo \\
& =  \lp\frac{b(n/k)}{R_{(k)}} \rp^2 \lmo \frac{1}{k}\sum_{i=1}^n \lp \frac{R_i}{b(n/k)} \rp^{2}   I_{R_i \ge R_{(k)} } - \frac{1}{k}\sum_{i=1}^n \lp \frac{R_i}{b(n/k)} \rp^{2}   I_{R_i \ge b(n/k) } \rmo \\
& =  \lmo \frac{1}{k}\sum_{i=1}^n \lp \frac{R_i}{R_{(k)}} \rp^{2}   I_{R_i \ge R_{(k)} } - \lp\frac{b(n/k)}{R_{(k)}} \rp^2  \frac{1}{k}\sum_{i=1}^n \lp \frac{R_i}{b(n/k)} \rp^{2}   I_{R_i \ge b(n/k) } \rmo
\end{align*}

Then, by Lemma <ref> (iii), we have that $\lp {b(n/k)}/{R_{(k)}} \rp^2 \convP 1$. By Lemma <ref> that $\frac{1}{k}\sum_{i=1}^n \lp{R_i}/{R_{(k)}} \rp^{2}   I_{R_i \ge R_{(k)} } \convP \ag/(\ag-2)$ and $\frac{1}{k}\sum_{i=1}^n \lp {R_i}/{b(n/k)} \rp^{2}   I_{R_i \ge b(n/k) }\convP \ag/(\ag-2)$. Therefore, we have that  $P_1(n) \convP 0$.

Now we work on $P_2(n)$. Observe that

\begin{align*}
P_2(n) & = \lmo \frac{1}{k}\sum_{i=1}^n \lip \frac{X_i}{R_i}, \frac{Y_i}{R_i} \rip R_i^2 \lp \frac{1}{R_{(k)}^2} - \frac{1}{b(n/k)^2}\rp
I_{R_i \ge b(n/k) } \rmo \\
& \le  \lmo \frac{b(n/k)^2}{R_{(k)}^2} - 1 \rmo  \frac{1}{k}\sum_{i=1}^n  \lmo \lip \frac{X_i}{R_i}, \frac{Y_i}{R_i}  \rip \rmo \lp \frac{R_i}{b(n/k)} \rp^{2}
I_{R_i \ge b(n/k)} \\
& \le  \lmo \frac{b(n/k)^2}{R_{(k)}^2} - 1 \rmo  \frac{1}{k}\sum_{i=1}^n \lp \frac{R_i}{b(n/k)} \rp^{2} I_{R_i \ge b(n/k)}.
\end{align*}

By Lemma <ref>, we have that $\frac{1}{k}\sum_{i=1}^n \lp {R_i}/{b(n/k)} \rp^{2} I_{R_i \ge b(n/k)} = O_P(1)$, and by Lemma <ref> (iii), we have that $b(n/k) / R_{(k)} \convP 1$. Thus, $P_2(n) \convP 0$.

Proof of Theorem <ref>. It follows from Propositions <ref> and <ref>.

§ PROOF OF LEMMA <REF> IN SECTION <REF>

We begin by noting that since $Z_1$ and $Z_2$ are independent, there exists $\nu$ in $M_+(\mbR_+^2)$ such that

\begin{equation} \label{e:dRV2}
nP \lp \frac{(|Z_1|, |Z_2|)}{b(n)} \in \cdot \rp \convv \nu,
\end{equation}

and for $\bx = [x_1, x_2]^{\top}$

\[
\nu([0, \bx]^c) = c\{(x_1)^{-\ag} +(x_2)^{-\ag}\}.
\]

With the choice of $b(n)$ defined by

\begin{align} \label{e:bn}
n^{-1} &= P(\| Z_1 \phi_1 \|  \vee  \| \rho Z_1 \phi_1+\sqrt{1-\rho^2} Z_2 \phi_2 \| > b(n))\\ \nonumber
&= P(|Z_1|  \vee  (\rho^2 Z_1^2 + (1-\rho^2) Z_2^2)^{1/2} > b(n)),
\end{align}

we set $c= 1/(1+ (1-\rho^2)^{\ag/2})$ to ensure that $\nu$ is a probability measure on $\{(z_1,z_2): |z_1| \vee ( \rho^2 z_1^2 + (1-\rho^2) z_2^2 )^{1/2} >1 \}$.

We claim that

\begin{align} \label{e:sxy}
& \sg_{XY}=\rho  \frac{c\ag}{\ag-2};\\ \label{e:sx}
& \sg_{X}^2 =\frac{c\ag}{\ag-2};\\ \label{e:sy}
& \sg_{Y}^2 =\lbr \rho^2 + (1-\rho^2)^{\ag/2} \rbr \frac{c\ag}{\ag-2}.
\end{align}

We first work on sxy. Since the terms with the $N_j$ do not affect the extremal behavior of $X$ and $Y$, we have that by Proposition <ref>

\begin{align*}
&\sg_{XY} \\
&= \lim_{n \to \infty} E \lb \lip \frac{Z_1 \phi_1}{b(n)}, \frac{\rho Z_1 \phi_1+\sqrt{1-\rho^2} Z_2 \phi_2}{b(n)}\rip  \Bigg| \| Z_1 \phi_1 \|  \vee  \| \rho Z_1 \phi_1+\sqrt{1-\rho^2} Z_2 \phi_2 \| > b(n)\rb  \\
&= \lim_{n \to \infty} \frac{1}{P(\| Z_1 \phi_1 \|  \vee  \| \rho Z_1 \phi_1+\sqrt{1-\rho^2} Z_2 \phi_2 \| > b(n))} \times \\
&\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ E \lb \lip \frac{Z_1 \phi_1}{b(n)}, \frac{\rho Z_1 \phi_1+\sqrt{1-\rho^2} Z_2 \phi_2}{b(n)}\rip  I_{\| Z_1 \phi_1 \|  \vee  \| \rho Z_1 \phi_1+\sqrt{1-\rho^2} Z_2 \phi_2 \| > b(n)}\rb  \\
&= \lim_{n \to \infty} \frac{1}{P(\| Z_1 \phi_1 \|  \vee  \| \rho Z_1 \phi_1+\sqrt{1-\rho^2} Z_2 \phi_2 \| > b(n))}   E \lb  \rho \frac{Z_1^2 }{b(n)^2} I_{|Z_1| \vee (\rho^2 Z_1^2 +(1-\rho^2)Z_2^2 )^{1/2} > b(n)} \rb
\end{align*}

It then follows from dRV2 and bn  that

\begin{align*}
\sg_{XY} &= \lim_{n \to \infty} n E \lb  \rho \frac{Z_1^2 }{b(n)^2} I_{|Z_1| \vee (\rho^2 Z_1^2 +(1-\rho^2)Z_2^2 )^{1/2} > b(n)} \rb \\
&= \lim_{n \to \infty} \int_{\mbR^2_+} \rho z_1^2 I_{ |z_1| \vee ( \rho^2 z_1^2 + (1-\rho^2) z_2^2 )^{1/2} >1} nP\lp \frac{|Z_1|}{b(n)} \in dz_1, \frac{|Z_2|}{b(n)}\in dz_2 \rp \\
&= \int_{\mbR_+^2}  \rho z_1^2 I_{ |z_1| \vee ( \rho^2 z_1^2 + (1-\rho^2) z_2^2 )^{1/2} >1} \ \nu(dz_1, dz_2) \\
& = \int_{\mbR_+}  \rho z_1^2 I_{ \{(z_1, 0): z_1>1\} } \ c\nu_{\ag}(dz_1) + \int_{\mbR_+}  \rho z_1^2 I_{ \{(0, z_2): z_2>1/(1-\rho^2)^{1/2} \}} \ c\nu_{\ag}(dz_2) \\
&= \int_{1}^{\infty}  \rho z_1^2 \ c\nu_{\alpha}(dz_1) + 0 = \rho \frac{c\ag}{\ag-2}.
\end{align*}

Analogously, for sx we can show that

\begin{align*}
&\sg_{X}^2 \\
&= \lim_{n \to \infty} E \lb \lip \frac{Z_1 \phi_1}{b(n)}, \frac{Z_1 \phi_1}{b(n)}\rip  \Bigg| \| Z_1 \phi_1 \|  \vee  \| \rho Z_1 \phi_1+\sqrt{1-\rho^2} Z_2 \phi_2 \| > b(n)\rb  \\
&= \lim_{n \to \infty} \frac{1}{P(\| Z_1 \phi_1 \|  \vee  \| \rho Z_1 \phi_1+\sqrt{1-\rho^2} Z_2 \phi_2 \| > b(n))}  E \lb  \frac{Z_1^2 }{b(n)^2} I_{|Z_1| \vee (\rho^2 Z_1^2 +(1-\rho^2)Z_2^2 )^{1/2} > b(n)} \rb \\
&= \lim_{n \to \infty} nE \lb  \frac{Z_1^2 }{b(n)^2} I_{|Z_1| \vee (\rho^2 Z_1^2 +(1-\rho^2)Z_2^2 )^{1/2} > b(n)} \rb \\
\end{align*}

Next, we work on sy. Observe that

\begin{align*}
&\sg_{Y}^2 \\
&= \lim_{n \to \infty} E \lb   \frac{\|\rho Z_1 \phi_1+\sqrt{1-\rho^2} Z_2 \phi_2 \|^2}{b(n)^2}  \Bigg| \| Z_1 \phi_1 \|  \vee  \| \rho Z_1 \phi_1+\sqrt{1-\rho^2} Z_2 \phi_2 \| > b(n)\rb  \\
&= \lim_{n \to \infty} \frac{1}{P(\| Z_1 \phi_1 \|  \vee  \| \rho Z_1 \phi_1+\sqrt{1-\rho^2} Z_2 \phi_2 \| > b(n))} \times \\
&\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ E \lb  \frac{\rho^2 Z_1^2 +(1-\rho^2) Z_2^2 }{b(n)^2} I_{|Z_1| \vee (\rho^2 Z_1^2 +(1-\rho^2)Z_2^2 )^{1/2} > b(n)} \rb.
\end{align*}

Then, again it follows from dRV2 and bn    that

\begin{align*}
\sg_{Y}^2 &= \lim_{n \to \infty} nE \lb  \frac{\rho^2 Z_1^2 +(1-\rho^2) Z_2^2 }{b(n)^2} I_{|Z_1| \vee (\rho^2 Z_1^2 +(1-\rho^2)Z_2^2 )^{1/2} > b(n)} \rb   \\
&= \lim_{n \to \infty} \int_{\mbR^2_+} \lbr \rho^2 z_1^2 + (1-\rho^2)z_2^2 \rbr I_{ |z_1| \vee ( \rho^2 z_1^2 + (1-\rho^2) z_2^2 )^{1/2} >1} nP\lp \frac{|Z_1|}{b(n)} \in dz_1, \frac{|Z_2|}{b(n)}\in dz_2 \rp \\
&= \int_{\mbR_+^2}  \lbr \rho^2 z_1^2 + (1-\rho^2)z_2^2 \rbr I_{ |z_1| \vee ( \rho^2 z_1^2 + (1-\rho^2) z_2^2 )^{1/2} >1} \ \nu(dz_1, dz_2) \\
& = \int_{\mbR_+}  \rho^2 z_1^2 I_{ \{(z_1, 0): z_1>1\} } \ c\nu_{\ag}(dz_1) + \int_{\mbR_+}  (1-\rho^2)z_2^2 I_{ \{(0, z_2): z_2>1/(1-\rho^2)^{1/2} \}} \ c\nu_{\ag}(dz_2) \\
&= \int_{1}^{\infty}  \rho^2 z_1^2 \ c\nu_{\alpha}(dz_1) + \int_{1/(1-\rho^2)^{1/2}}^{\infty}   (1-\rho^2)z_2^2 \ c\nu_{\alpha}(dz_1) \\
&= \lbr \rho^2  + (1-\rho^2)^{\ag/2} \rbr \frac{c\ag}{\ag-2}.
\end{align*}

§ CONSISTENCY OF $\HAT{\RHO}_{N,K}$ FOR $\AG \IN \{3,4,5\}$

Empirical biases (standard errors) of $\hat{\rho}_{n,k}$ when $\ag=3$. The KS method is used to choose optimal $k$s. On average, it selects $k=8 \sim 9$ ($N=100$) and $k=29 \sim 33$ ($N=500$).

$\rho_{XY}$     $N=100$      $N=500$     $\rho_{XY}$     $N=100$      $N=500$

0        -0.002 (0.117)       0.002 (0.062)

0.1      -0.011 (0.133)       0.000 (0.081)        -0.1      0.002 (0.130)      0.003 (0.076)

0.2      -0.009 (0.167)     -0.012 (0.109)        -0.2      0.014 (0.156)     0.008 (0.116)

0.3      -0.045 (0.185)      -0.020 (0.130)       -0.3      0.035 (0.187)     0.022 (0.124)

0.4      -0.061 (0.193)      -0.041 (0.134)       -0.4      0.055 (0.194)     0.036 (0.145)

0.5      -0.097 (0.202)      -0.052 (0.150)       -0.5      0.084 (0.190)     0.064 (0.149)

0.6      -0.127 (0.201)      -0.078 (0.159)       -0.6      0.126 (0.211)     0.080 (0.152)

0.7      -0.156 (0.197)      -0.106 (0.148)       -0.7      0.161 (0.196)     0.108 (0.150)

0.8      -0.183 (0.184)      -0.130 (0.135)       -0.8      0.194 (0.195)     0.118 (0.141)

0.9      -0.219 (0.165)      -0.142 (0.122)       -0.9      0.222 (0.168)     0.142 (0.121)

1.0      -0.222 (0.125)      -0.139 (0.079)       -1.0      0.230 (0.132)     0.141 (0.082)

Empirical biases (standard errors) of $\hat{\rho}_{n,k}$ when $\ag=4$. The KS method is used to choose optimal $k$s. On average, it selects $k=8$ ($N=100$) and $k=25 \sim 29$ ($N=500$).

$\rho_{XY}$     $N=100$      $N=500$     $\rho_{XY}$     $N=100$      $N=500$

0       0.000 (0.172)        0.003 (0.094)

0.1      -0.030 (0.167)       -0.019 (0.093)       -0.1      0.027 (0.163)     0.019 (0.097)

0.2      -0.053 (0.188)      -0.052 (0.114)        -0.2      0.064 (0.184)     0.051 (0.116)

0.3      -0.123 (0.198)      -0.085 (0.127)        -0.3      0.102 (0.180)     0.083 (0.130)

0.4      -0.157 (0.192)        -0.118 (0.137)       -0.4      0.147 (0.202)     0.121 (0.136)

0.5      -0.211 (0.206)       -0.155 (0.150)       -0.5      0.197 (0.196)     0.163 (0.140)

0.6      -0.253 (0.202)       -0.204 (0.146)       -0.6      0.256 (0.199)     0.209 (0.154)

0.7      -0.306 (0.202)      -0.249 (0.156)        -0.7      0.315 (0.205)     0.255 (0.154)

0.8      -0.366 (0.207)       -0.291 (0.150)       -0.8      0.366 (0.206)     0.291 (0.150)

0.9      -0.415 (0.193)       -0.340 (0.141)       -0.9      0.415 (0.200)     0.329 (0.142)

1.0      -0.410 (0.170)       -0.326 (0.124)       -1.0      0.428 (0.171)     0.325 (0.126)

Empirical biases (standard errors) of $\hat{\rho}_{n,k}$ when $\ag=5$. The KS method is used to choose optimal $k$s. On average, it selects $k=7 \sim 8$ ($N=100$) and $k=18 \sim 20$ ($N=500$).

$\rho_{XY}$     $N=100$      $N=500$     $\rho_{XY}$     $N=100$      $N=500$

0       -0.007 (0.199)       -0.004 (0.159)

0.1      -0.046 (0.208)       -0.043 (0.144)       -0.1      0.053 (0.195)     0.041 (0.145)

0.2      -0.087 (0.204)       -0.095 (0.145)       -0.2      0.097 (0.202)     0.096 (0.158)

0.3      -0.173 (0.210)      -0.147 (0.149)       -0.3      0.155 (0.201)     0.142 (0.161)

0.4      -0.229 (0.198)       -0.203 (0.151)       -0.4      0.212 (0.209)     0.196 (0.161)

0.5      -0.289 (0.202)       -0.263 (0.164)        -0.5      0.279 (0.206)     0.264 (0.156)

0.6      -0.351 (0.217)       -0.324 (0.172)       -0.6      0.356 (0.207)     0.330 (0.184)

0.7      -0.411 (0.206)       -0.388 (0.171)       -0.7      0.424 (0.210)     0.397 (0.157)

0.8      -0.491 (0.209)      -0.458 (0.170)        -0.8      0.491 (0.207)     0.446 (0.162)

0.9      -0.556 (0.207)       -0.511 (0.156)       -0.9      0.554 (0.207)     0.505 (0.163)

1.0      -0.548 (0.171)       -0.502 (0.143)       -1.0      0.564 (0.178)     0.509 (0.149)