# Treatment Effects in Extreme Regimes

Ahmed Aloui1 Ali Hasan2 Yuting Ng1 Miroslav Pajic1 Vahid Tarokh1
1 Department of Electrical and Computer Engineering  Duke University
2 Department of Biomedical Engineering  Duke University
<EMAIL_ADDRESS>
###### Abstract

Understanding treatment effects in extreme regimes is important for
characterizing risks associated with different interventions. This is hindered
by the unavailability of counterfactual outcomes and the rarity and difficulty
of collecting extreme data in practice. To address this issue, we propose a
new framework based on extreme value theory for estimating treatment effects
in extreme regimes. We quantify these effects using variations in tail decay
rates of potential outcomes in the presence and absence of treatments. We
establish algorithms for calculating these quantities and develop related
theoretical results. We demonstrate the efficacy of our approach on various
standard synthetic and semi-synthetic datasets.

##  1 Introduction

Understanding the effect of an intervention is a primary objective of causal
inference. In the Neyman-Rubin framework [Neyman, 1923, Rubin, 1974], the
average treatment effect (ATE) and the conditional average treatment effect
(CATE) are the quantities of general interest. These metrics are instrumental
for informed decision-making [Abadie and Cattaneo, 2018] and have been widely
applied to various areas such as econometrics and epidemiology [Rothman and
Greenland, 2005].

Despite their importance, these statistics are based on the _expectation_ of
the difference of the outcome variables influenced by the intervention.
However, since they are focused on the expectation, these quantities do not
consider the extreme behavior of the treatment effect. In many scenarios, it
is equally or perhaps more important to understand the _worst case_ behavior
of interventions rather than only the expected behavior. For example, we may
prescribe a policy that is beneficial on average, yet it increases the
likelihood of extremely adverse events that we want to completely prevent.
This is particularly relevant at an individual level, where an intervention
can lead to severe adverse effects for a certain subset of members of the
population compared to others. Applied in a medical setting, a specific
patient may experience stronger adverse side effects as a response to
treatment than those observed on the population level or in the average case.
In turn, there is a need to characterize the _maximum_ (or similarly the
_minimum_) of the treatment effect rather than the average.

This provides the motivation for the present work where we develop a framework
for estimating treatment effects in the _tails_ of the distribution. Although
quantile methods have found application in causal inference [Díaz, 2017,
Giessing and Wang, 2021], their performance becomes notably unstable and
inaccurate for high quantile statistics, especially when data is limited.
Therefore, we use the extrapolation properties of Extreme Value Theory (EVT)
[Haan and Ferreira, 2006] to propose a measure quantifying extreme treatment
effects. Under certain assumptions, EVT provides asymptotic theory for the
behavior of data within the tail of a distribution through the Generalized
Extreme Value (GEV) distributions. GEVs have location and scale parameters as
well as a shape parameter, which is particularly important since it measures
the tail decay rate. We propose defining the extreme treatment effect as the
difference between the shape parameters of potential outcomes. This
intuitively quantifies the change in the distribution tail decay rate upon
treatment administration. We refer to these statistics as the _extreme
treatment effect_ (ETE). Related to this, we also examine the differences in
the tails of the individual level potential outcomes referred to as the
_conditional extreme treatment effect_ (CETE).

##### Motivating Example

In applications such as healthcare, it may be necessary to ensure that the
tail behavior of an outcome (e.g. corresponding to adverse reactions to a
drug) does not significantly change under the treatment intervention. We can
use the ETE to characterize which distribution is likely to have more
realizations deep in the tail since this quantity describes the difference in
the decay rate of the tails. Faster tail decay will correspond to fewer
extreme events and vice-versa. The idea is illustrated in Figure 1 where the
outcome above $12$ corresponds to a serious adverse event but below is not
considered harmful. On average, the treatment group appears more beneficial
(since it has a negative ATE), but it simultaneously induces a greater risk of
serious adverse effects (and thus has a higher ETE). Additionally, we focus on
estimating the CETE for the _individual_ level rather than the population
level, which has implications in fields such as personalized medicine. We
illustrate this idea in greater detail in Appendix A. Note that for the
problems that we are studying, higher order moments such as variance or skew
are insufficient for understanding the tail risk. The reason is due to the
fact that many of the underlying phenomena are heavy tailed, and higher order
moments may not be defined for these data, resulting in statistics that never
converge.

Figure 1: Illustration of what the ETE represents. The control group has
faster decaying tails with a shape parameter of $0.5$ whereas the treatment
group has a shape parameter of $1.5$ with slowly decaying tails. Yet, the
treatment group has a lower mean. Hence, ETE $=1$, while the ATE is negative.

While capturing essential characteristics of the potential outcomes, the
proposed ETE and CETE have challenges associated with estimation from data.
The first challenge is due to the lack of observed data in the tails of the
distributions. By definition, data in the tails are rare, hence we would not
expect to have access to many samples from the tail. To counteract this
challenge, our framework uses EVT to represent the tail distribution as a
parametric class of distributions and extrapolate deep into the tails from
near-tail data. The second challenge arises from the unavailability of
counterfactual outcomes, which is referred to as the _fundamental problem of
causal inference_ [Rubin, 1974, Holland, 1986]. This challenge persists in any
causal inference setting and is due to the fact that we cannot observe the
outcome for the same individual under two different treatment regimes. The
third challenge is due to _unobserved confounders_. While it is common to
assume that the observations are generated such that no unobserved confounders
exist, which is known as the unconfoundedness assumption111Also referred to as
conditional exchangeability., it is often very hard to test this assumption in
practice.

To address these challenges, we propose a new _tail unconfoundedness_
condition that we prove to be strictly weaker than unconfoundedness. Moreover,
we develop statistical estimators for the proposed ETE and CETE and we prove
identifiability and consistency. We illustrate these results numerically in
various experiments with both synthetic and semi-synthetic datasets.

Related Work As discussed in the introduction, EVT provides a mathematical
theory for extrapolating to the tails of a distribution [Haan and Ferreira,
2006]. In causal inference, methods for estimating the Average Treatment
Effect (ATE) and Conditional Average Treatment Effect (CATE), extreme average
treatment effects, and quantile treatment effects have been explored in the
literature. Methods for estimating the ATE include covariate adjustment of
Rubin [1974], also known as the backdoor adjustment in the do-calculus
framework [Pearl, 2009]. With regards to deep learning methods for estimating
CATE, the $s$-learner and $t$-learner neural networks [Künzel et al., 2019],
counterfactual regression [Johansson et al., 2016], and TARNet [Shalit et al.,
2017] are well-known examples. Additional machine learning approaches such as
random forests [Wager and Athey, 2018] and Gaussian processes [Alaa and Van
Der Schaar, 2017] have also been employed for CATE estimation. Other methods
have developed doubly robust learners for CATE estimation [Kennedy, 2020].
Estimators for extreme quantile treatment effects were developed in [Zhang,
2018], but without considering extrapolation outside the range of
observations. Inferring the conditional value at risk of a treatment effect
was investigated in Kallus [2023] but the case of tail quantiles was not
explicitly considered. Estimating the extreme average treatment effect was
addressed in [Huang et al., 2022] using extreme value theory. Works on
estimating a parametric or semi-parametric form of the counterfactual
distribution include Kennedy et al. [2021], Chernozhukov and Hansen [2005],
Chernozhukov et al. [2013]. Methods for estimating the quantiles include
[Díaz, 2017, Firpo, 2007, Frölich and Melly, 2013, Zhang et al., 2012,
Giessing and Wang, 2021]. These methods do not consider estimating the
conditional extreme treatment effect and focus on the population level.
Moreover, we provide a different framework to define these quantities of
interest and provide a detailed mathematical motivation that exploits the
extrapolating behavior of EVT. Our approach is closely related to the study of
max-stable processes, which have been applied across various disciplines to
model extremes over function spaces [Davison et al., 2012]. The key difference
is these distributions are often estimated assuming an abundance of temporal
observations, which is an assumption that may not hold in causal inference
settings.

##  2 Background

###  2.1 Causal inference

Consider covariates $X\in\mathcal{X}$ with $\mathcal{X}\subset\mathbb{R}^{d}$
being the domain. Let $T\in\mathcal{T}=\\{0,1\\}$ denote the treatment
variable. The group assigned $T=0$ and $T=1$ are called the control and
treatment groups, respectively. Denote the variable $Y\in\mathbb{R}$ as the
factual outcome. Potential outcomes are the outcomes that would have been
observed if only treatment $T=1$ or $T=0$ was assigned and are respectively
represented by $Y_{1}$ and $Y_{0}$ [Morgan and Winship, 2015]. The
relationship between the factual outcome and the potential outcomes is
$Y=Y_{1}T+Y_{0}(1-T)$. The main challenge in causal inference is that we do
not observe the counterfactual outcomes, i.e., the outcomes that would have
been observed if the treatment were reversed. We call any statistics of the
potential outcomes a _causal statistics_. A causal statistic is said to be
identifiable if it can be reduced to a function of the factual outcomes $Y$,
the treatment assignment $T$, and the covariates $X$ [Imbens, 2004]. This is
formalized in the following definition:

######  Definition 2.1 (Causal Identifiability).

A causal statistic (e.g $\mathbb{E}[Y_{1}]$) is identifiable if it can be
computed from a purely observational quantity (e.g.
$\mathbb{E}[\mathbb{E}[Y\mid T=1,X]]$).

Similarly, an estimator, as a particular statistic, is said to be identifiable
when it meets this criterion. We next describe some assumptions that are often
made in order to guarantee the causal identifiability of CATE and ATE. The
first is consistency, which requests that the potential outcomes align with
the factual outcome:

######  Assumption 2.2 (Consistency).

For $t\in\mathcal{T}$, the following holds, $T=t\implies Y_{t}=Y$.

The next is overlap222Also known as the positivity assumption. which requires
that the different treatment groups have the same support:

######  Assumption 2.3 (Overlap).

$\forall x\in\mathcal{X},\;0<P(T=1\mid X=x)<1$

Finally we assume unconfoundedness, which requires that the potential outcomes
are independent of the treatment given the covariates:

######  Assumption 2.4 (Unconfoundedness).

$\left(Y_{1},Y_{0}\right)\perp\\!\\!\\!\\!\perp T\mid X$

We will later relax this assumption to a weaker _tail unconfoundedness_
assumption, which intuitively states that this property only needs to hold in
the tail of the outcome variable.

###  2.2 Extreme value theory

EVT provides a principled mathematical framework for extrapolating to the
tails of a distribution. The Fisher-Tippett-Gnedenko theorem states that if
the shifted and scaled versions of maximum order statistics of a distribution
$P(\cdot)$ converge to a non-degenerate distribution, then it is guaranteed to
converge to a family of well-specified distributions known as the generalized
extreme value (GEV) distributions [Haan and Ferreira, 2006, Fisher and
Tippett, 1928, Gnedenko, 1943]. If this convergence holds, the distribution
$P(\cdot)$ is said to be in the _maximum domain of attraction (MDA)_ 333The
MDA of a distribution $p$ is the distribution to which the maximum of samples
from $p$ converges to. of the GEV. The distribution is characterized by three
parameters: $\mu$, the location parameter; $\sigma$, the scale parameter; and,
$\xi$ the _shape_ or _tail index_ parameter.

######  Definition 2.5.

The GEV distribution is defined as:

$G_{\xi}(y)=\exp\left(-(1+\xi\bar{y})^{-1/\xi}\right),\>\>1+\xi\bar{y}>0$ (1)

where $\bar{y}=\frac{y-\mu}{\sigma}$ with parameters $\mu,\sigma,\xi$ The case
$\xi=0$ is defined as the limit where $\xi$ approaches zero.

For our purposes the shape parameter $\xi$ is of interest, which governs the
rate of tail decay. When $\xi>0$ the distribution exhibits a heavy tail,
whereas when $\xi<0$ the tail is bounded. If $\xi=0$, the tail decays
exponentially. The Fréchet, Weibull, and Gumbel distributions respectively
correspond to positive, negative, and zero values of $\xi$.

######  Theorem 2.6 (Fisher–Tippett–Gnedenko).

Let $Y^{(1)},Y^{(2)},\ldots,Y^{(n)}$ be a sequence of iid real random
variables. If there exist two sequences of real numbers $a_{n}>0$ and
$b_{n}\in\mathbb{R}$ such that the following limits converge to a non-
degenerate distribution function:

$\lim_{n\rightarrow\infty}P\left(\frac{\max\left\\{Y^{(1)},\ldots,Y^{(n)}\right\\}-b_{n}}{a_{n}}\leq
y\right)=G(y),$

then the limiting distribution $G$ is the GEV distribution.

In general, fitting such distributions requires many repeated observations
such that the maximum can be computed, thus obtaining samples approaching the
tail. The _block maxima_ method refers to the approach when the data is
subdivided into blocks and the maximum quantity per block is used for fitting
[Coles et al., 2001]. Specifically, $K\times n$ observations
$\left\\{Y^{(i,j)}\right\\}_{i,j=1}^{K\times n}$ observations are divided into
$K$ blocks of size $n$. Let $M^{(i)}=\underset{j=1,\ldots,K}{\max}Y^{(i,j)}$
for $i=1,\ldots,n$. The values $\left\\{M^{(i)}\right\\}_{i=1}^{n}$ are then
used to fit the GEV model. In the causal inference setting, multiple
observations per individual is often prohibitively expensive to obtain, so we
introduce alternatives to the block maxima method later in the text.

##  3 Theoretical Framework

We now present the proposed framework for estimating treatment effects in
extreme regimes. We begin by defining the estimation problem motivated by the
previous examples. Then, we establish an identifiability theorem under the
condition of tail unconfoundedness, which is a weaker requirement than
conditional unconfoundedness. Combining these results leads to a numerical
algorithm that can be readily implemented. The proofs of our results are
presented in Appendix B.

###  3.1 Problem setup

Given a dataset $S=\\{(x^{(i)},t^{(i)},z^{(i)})\\}_{i=1}^{N}$, sampled from a
$(X,T,Z)$ where $Z$ is the limiting tail distribution for $Y$. Our goal is to
estimate the effect of the treatment assignment on the tail indices of the
potential outcomes $Y_{0}$ and $Y_{1}$ as well as the conditioned potential
outcomes $Y_{0}|X=x$ and $Y_{1}|X=x$. This estimation faces the major
challenge in causal inference that the potential outcomes are only partially
observed and are subject to selection bias. To address this, we investigate
the identifiability of the conditional maximum likelihood estimator under tail
unconfoundedness. When presenting the method, we assume that the domain of
attraction of each individual can be estimated, and we assume repeated
observations for each individual. In practice, only one factual outcome per
individual is usually observed, making it difficult to estimate the maxima
through standard preprocessing methods such as the block maxima method. We
provide a practical way of dealing with this in Appendix C. Let
$x\in\mathcal{X}$ and consider the potential outcomes $Y_{1}$ and $Y_{0}$. The
outcomes conditioned on the individual $X=x$ are denoted as
$Y_{1}(x)=(Y_{1}\mid X=x)$ and $Y_{0}(x)=(Y_{0}\mid X=x)$. We assume both the
potential outcomes and their conditionals converge to non-degenerate
distributions within the domain of attraction of the Generalized Extreme Value
distributions:
$\text{GEV}(\stackunder[1.2pt]{$\mu$}{\rule{3.44444pt}{0.32289pt}}_{1},\stackunder[1.2pt]{$\sigma$}{\rule{3.44444pt}{0.32289pt}}_{1},\stackunder[1.2pt]{$\xi$}{\rule{3.44444pt}{0.32289pt}}_{1})$,
$\text{GEV}(\stackunder[1.2pt]{$\mu$}{\rule{3.44444pt}{0.32289pt}}_{0},\stackunder[1.2pt]{$\sigma$}{\rule{3.44444pt}{0.32289pt}}_{0},\stackunder[1.2pt]{$\xi$}{\rule{3.44444pt}{0.32289pt}}_{0})$,
$\text{GEV}(\mu_{1}(x),\sigma_{1}(x),\xi_{1}(x))$, and
$\text{GEV}(\mu_{0}(x),\sigma_{0}(x),\xi_{0}(x))$. Thus, the conditions of the
Fisher-Tippet-Gnedenko theorem are satisfied.

######  Definition 3.1 (ETE and CETE).

The extreme treatment effect is defined as:

$\stackunder[1.2pt]{\tau}{\rule{3.44444pt}{0.32289pt}}_{ext}=\stackunder[1.2pt]{\xi}{\rule{3.44444pt}{0.32289pt}}_{1}-\stackunder[1.2pt]{\xi}{\rule{3.44444pt}{0.32289pt}}_{0}.$

The conditional extreme treatment effect is defined as:

$\tau_{ext}(x)=\xi_{1}(x)-\xi_{0}(x).$

We introduce two metrics: $\varepsilon_{\text{ETE}}$ which quantifies the
error in estimating ETE, and $\varepsilon_{\text{CETE}}$, measuring the
precision in estimating CETE.

######  Definition 3.2.

Given estimators
$\hat{\stackunder[1.2pt]{$\tau$}{\rule{3.44444pt}{0.32289pt}}}_{ext}$ for ETE
and $\hat{\tau}_{ext}(x)$ for CETE, the performance errors are respectively
defined as:

$\varepsilon_{\text{ETE}}=|\hat{\stackunder[1.2pt]{\tau}{\rule{3.44444pt}{0.32289pt}}}_{ext}-\stackunder[1.2pt]{\tau}{\rule{3.44444pt}{0.32289pt}}_{ext}|$

and

$\varepsilon_{\text{CETE}}=\mathbb{E}\left[\left(\hat{\tau}_{ext}(x)-\tau_{ext}(x)\right)^{2}\right].$

For our proposed methods, we will use this precision measure for quantifying
the performance of different methods in our empirical studies.

###  3.2 Mathematical Motivation

In this section, we provide two mathematical interpretations of the ETE and
CETE to motivate why these quantities are interesting in a causal inference
setting.

##### Regularly Varying Perspective.

In this perspective, we can relate the ETE and CETE to the different rates of
tail decays of a function through regular variation. The basic interpretation
is that when both tail indices are positive, the difference between them
relates to the ratio of the survival distributions and can be used to study
which group is more likely to have outcomes exceed a certain threshold. We
review some definitions on regularly varying functions from Resnick [2007,
Chapter 1] and refer the reader to the chapter for additional information. A
function $U:\mathbb{R}_{+}\rightarrow\mathbb{R}_{+}$ is regularly varying
($RV$) if $\exists\,\xi\in\mathbb{R}$ such that

$\forall x>0,\lim_{t\rightarrow\infty}\frac{U(tx)}{U(t)}=x^{\xi}.$

This definition is due to Karamata [Bingham et al., 1989]. The parameter $\xi$
is called the exponent of variation, and we write $U\in RV(\xi)$ to mean that
$U$ is a regularly varying function with exponent of variation $\xi$. We will
present the results for the Fréchet case where $\xi>0$. The following remark
characterizes our motivation Resnick [1988]

######  Remark 3.3.

Let
$\stackunder[1.2pt]{$\xi$}{\rule{3.44444pt}{0.32289pt}}_{0},\stackunder[1.2pt]{$\xi$}{\rule{3.44444pt}{0.32289pt}}_{1}>0$.
Let $S_{1}(Y_{1}),S_{0}(Y_{0})$ be the survival distributions of the potential
outcome variables with positive support for the treatment and control groups
respectively. Then if
$\stackunder[1.2pt]{$\tau$}{\rule{3.44444pt}{0.32289pt}}_{ext}>0$, then
${S_{0}(Y_{0})}/{S_{1}(Y_{1})}\in
RV(\stackunder[1.2pt]{$\xi$}{\rule{3.44444pt}{0.32289pt}}_{1}-\stackunder[1.2pt]{$\xi$}{\rule{3.44444pt}{0.32289pt}}_{0})$.
Moreover, $S_{1}(Y_{1})<S_{0}(Y_{0})$ for $Y_{t}\gg 1$.

The proof for this remark relies on the Karamata representation theorem
[Galambos and Seneta, 1973], with further details in Appendix B. Therefore,
the ETE and CETE can be interpreted as characterizing the ratio of the
survival functions between the treatment and control groups.

##### Point Process Perspective.

Another practical viewpoint of EVT involves interpreting the GEV in terms of
the intensity of a Poisson point process [Coles et al., 2001, Chapter 7]. The
intensity of the Poisson process is given by the negative log of the density
in (1). One can use this form to calculate the expected number of extreme
points within a region of the state space. We formalize this in the following
remark:

######  Remark 3.4 (Expected Number of Events).

Consider a causal inference task with potential outcomes variable $Y_{t}$,
$t\in\\{0,1\\}$, that $Y_{t}$ is in the maximum domain of attraction of a GEV
with rate $\stackunder[1.2pt]{$\xi$}{\rule{3.44444pt}{0.32289pt}}_{t}$, and
ETE given by
$\stackunder[1.2pt]{$\tau$}{\rule{3.44444pt}{0.32289pt}}_{ext}>0$. Suppose
that $\mathcal{R}\subset\mathbb{R}$ is the Borel set of interest that should
be quantified for the causal inference task with $|\inf\mathcal{R}|\gg 0$.
Then, $\mathbb{E}[N_{Y_{1}}(\mathcal{R})]>\mathbb{E}[N_{Y_{0}}(\mathcal{R})]$
where
$N_{Y_{t}}(\mathcal{R})=\sum_{i=1}^{\infty}\delta_{Y_{t}^{(i)}\in\mathcal{R}}$
for $t\in\\{0,1\\}$ is a counting process over outcomes on some Borel subset
of $\mathbb{R}$.

The analog of this remark can easily be made for the CETE case by conditioning
on the individual covariates. This is particularly important, for example, in
a medical setting where one may want to consider the number of severe adverse
events in a treatment group versus a control group. If the outcome variable
being above a certain threshold corresponds to a serious adverse event, then
the ETE and CETE can be used to quantify the number of events.

###  3.3 Tail unconfoundedness and identifiablity

In order for the problem of estimating ETE and CETE to be tractable, we need
to establish identifiability of the statistics involving counterfactual
quantities in Definition 2.1 before these quantities can be estimated from
observational quantities. To this end, in addition to Assumptions 2.2 and 2.3,
we will introduce a new assumption in order to guarantee identifiability.

######  Assumption 3.5 (Tail Unconfoundedness).

$\forall t\in\mathcal{T},\;{P(Z_{t}\mid X,T)=P(Z_{t}\mid X)}$

where $Z_{t}$ is the MDA of the potential outcome $Y_{t}$.

Note that Assumption 3.5 is a weaker condition than the unconfoundedness
presented in Assumption 2.4, which is assumed in ATE and CATE estimation
[Imbens, 2004]. Assumption 2.4 requires that all the confounding variables are
observed. On the other hand, Assumption 3.5 only requires observing the
confounding variables that affect the tail. This is an important relaxation
since, even if some covariates that do not affect the tail are missing, ETE
and CETE can still be estimated.

##### Example where tail unconfoundedness is strictly weaker than
unconfoundedness.

To illustrate how the relaxation in the assumption is important, we present
discrete and continuous examples of distributions that satisfy Assumption 3.5
but not Assumption 2.4. For the discrete example let $P(Y_{1}\mid T=t,X)$ be a
Bernoulli distribution with a non-binary treatment parameter
$t\in\left(0,1\right)\\}$. For the continuous example let $P(Y_{1}\mid
X=x,T=t)$ follow a Beta distribution $\mathcal{B}(\alpha,\beta)$ with
parameters $\alpha=x$ and $\beta=t+1$ for $x,t>0$. Additional details are
provided in Appendix D. Moreover, Assumption 2.4 implies Assumption 3.5
formalized in the proposition below.

######  Proposition 3.6.

If unconfoundedness (Assumption 2.4) holds then so does the tail
unconfoundedness (Assumption 3.5).

Next, we define our conditional maximum likelihood estimator for estimating
the individual tail indices. We assume that the location and the scale
parameters are known. Let
$\mathcal{H}:\mathcal{X}\times\mathcal{T}\to\mathbb{R}$ be a class of
hypothesis functions estimating the shape parameter $\xi_{t}(x)$ for
$x,t\in\mathcal{X}\times\mathcal{T}$. In this section, we will assume that the
location and scale parameters $\mu,\sigma$ are known. In our empirical
results, $\mathcal{H}$ is considered to be a class of neural networks with
parameters $\theta\in\Theta$, i.e
$\mathcal{H}=\\{g_{\theta}:\mathcal{X}\times\mathcal{T}\to\mathbb{R}|\theta\in\Theta\\}$.

######  Definition 3.7 (CETE Estimator).

For $t\in\\{0,1\\}$, let $Z_{t}(x)$ be the conditional random variable
$Z_{t}|X=x$ for the maximum domain of attraction of $Y_{t}|X=x$. The shape
parameter $\xi_{t}(x)$ estimator is defined as:

$\theta^{*}=\arg\max_{\theta\in\Theta}\sum_{j=1}^{N}\log
p\left(Z_{t}(x^{(j)})=z_{t}^{(j)}\;\big{|}\;g_{\theta}\left(x^{(j)},t\right)\right)$

and for every $x\in\mathcal{X}$,

$\hat{\xi}_{t}(x)=g_{\theta^{*}}(x,t)$

Accordingly, we define our CETE estimator as:
$\hat{\tau}_{ext}(x)=\hat{\xi}_{1}(x)-\hat{\xi}_{0}(x).$

Since $Z_{1}|X=x$ is given as the tail distribution, $p$ is taken to be the
density of the conditional GEV in Definition 2.5. A naive estimator for the
ETE can be defined as follows:

$\tilde{\stackunder[1.2pt]{\tau}{\rule{3.44444pt}{0.32289pt}}}=\tilde{\stackunder[1.2pt]{\xi}{\rule{3.44444pt}{0.32289pt}}}_{1}-\tilde{\stackunder[1.2pt]{\xi}{\rule{3.44444pt}{0.32289pt}}}_{0}.$

where,

$\tilde{\stackunder[1.2pt]{\xi}{\rule{3.44444pt}{0.32289pt}}}_{t}=\arg\max_{\phi}\sum_{j=1}^{N}\log
p\left(Z=z^{(j)}|T=t;\phi\right)$

Therefore, the naive estimator takes the samples from the GEV of the outcomes
conditioned on the treatment and fits a GEV to it. This results in a biased
estimator of the ETE because of the selection bias introduced by the treatment
selection. Now we define our proposed ETE estimator based on the the CETE
estimator.

######  Definition 3.8 (ETE Estimator).

For $t\in\\{0,1\\}$, let $\\{\hat{\xi}_{t}(x^{(j)})\\}$ be the estimated
conditional shape parameter using the CETE estimator. We sample
$\\{\hat{z}^{(j)}_{t}\\}_{j=1}^{N}$ from the GEV with the estimated shape
parameters. The tail index
$\stackunder[1.2pt]{$\xi$}{\rule{3.44444pt}{0.32289pt}}_{t}$ estimator is
defined as:

$\widehat{\stackunder[1.2pt]{\xi}{\rule{3.44444pt}{0.32289pt}}}_{t}=\arg\max_{\phi}\sum_{j=1}^{N}\log
p\left(\hat{Z}_{t}=\hat{z}^{(j)}_{t};\phi\right)$

Accordingly, we define our ETE estimator as:
$\hat{\stackunder[1.2pt]{$\tau$}{\rule{3.44444pt}{0.32289pt}}}=\hat{\stackunder[1.2pt]{$\xi$}{\rule{3.44444pt}{0.32289pt}}}_{1}-\hat{\stackunder[1.2pt]{$\xi$}{\rule{3.44444pt}{0.32289pt}}}_{0}.$

The following theorem makes the defined causal quantities statistically
tractable.

######  Proposition 3.9 (Identifiability).

Under Assumptions 2.2, 2.3, and 3.5, the CETE and ETE estimators defined in
Definition 3.7 and Definition 3.8 respectively are causally identifiable.

The identifiability theorem reduces our original causal statistics (that
depend on tails of the potential outcomes) to quantities that only depend on
observable variables. However, even under an identifiability assumption, the
estimator in Definition 3.7 may not be easy to compute due to the fact that
only one observation per individual may exist. We include a theoretical
discussion as well as empirical results when only one outcome is observed per
individual in Appendix C.

We now prove the following statistical inference results for the proposed
estimators.

######  Remark 3.10.

The ETE estimator defined in Definitions 3.8 is consistent and asymptotically
normal.

This remark asserts that, with an increasing number of observed data points,
the ETE and CETE estimators converge to their true values and their
distributions approach a Gaussian.

##  4 Computational Method

In this section, we present the computational method employed in our study.
Given a dataset $\\{(x_{i},t_{i},y_{i})\\}_{i=1}^{N}$, we first apply a
variation of the block maxima method. We then learn the different parameters
of the limiting GEV by maximizing the conditional log-likelihood of the
maximum of the outcomes to the covariates. The CETE and ETE can then be
computed using the estimated quantities on $\xi$, a pseudocode of the
algorithm is provided in Appendix E.

Empirical $\varepsilon$-max-sampler. As noted, the standard procedure of using
the block maxima in time cannot be used for practical scenarios of interest in
the causal setting since we do not have a temporal component — we only observe
one outcome for every individual. To approach the task of approximating the
maximum of an outcome variable across multiple copies, we instead consider the
analogy of the block maxima in _space_ to exploit the property that
individuals with similar characteristics tend to have similar outcomes. As a
result, we can use realizations of individuals locally to approximate
independent and identical realizations of a single individual. As the number
of points increases and the ball size around every individual converges to
zero the approximation of multiple copies of individual outcomes becomes
increasingly accurate when taking the maximum of their outcomes.

To do this, we use the $k$-means algorithm to identify clusters of individuals
with similar feature profiles. Then, within each cluster, we select the
individual with the largest outcome as an approximation of the maximum outcome
for that group of similar individuals. We provide a detailed description of
the $\varepsilon$-max-sampler procedure in Appendix E as well as a theoretical
analysis of this empirical approach in Appendix C.

##### Likelihood-Based Estimator.

Since we assume that the potential outcomes are in the maximum domain of
attraction of a GEV, we maximize the likelihood of a GEV parameterized by
functions mapping the observations to the parameters. Let
$\theta_{t}=(\mu_{t},\sigma_{t},\xi_{t})$ where
$\mu_{t}:\mathcal{X}\to\mathbb{R}$, $\sigma_{t}:\mathcal{X}\to\mathbb{R}_{+}$
and $\xi_{t}:\mathcal{X}\to\mathbb{R}$ for $t\in\\{0,1\\}$. Let
$\\{(x^{(m_{i})},t^{(m_{i})},y^{(m_{i})})\\}_{i=1}^{K}$ $K$ be the data points
after computing the $\varepsilon$-max-sampler. Then we have for
$t\in\\{0,1\\}$,

$\displaystyle\hat{\theta}_{t}=\arg\max_{\theta}\sum_{i=1}^{K}\log(\ell(\theta\mid
x^{(m_{i})},t^{(m_{i})},y^{(m_{i})}))\mathbbm{1}_{(t^{m_{i}}=t)}$

where

$\ell(\theta\mid x,t,y)=\frac{\Phi(y\mid
x,t)^{\xi_{t}(x)+1}}{\sigma_{t}(x)}e^{-\Phi(y\mid x,t)}\;$

and

$\Phi(y\mid
x,t)=\left(1+\xi_{t}(x)\left(\frac{y-\mu_{t}(x)}{\sigma_{t}(x)}\right)\right)^{-1/\xi_{t}(x)}.$

After fitting the parameters, we can estimate the CETE as the following
$\hat{\tau}_{ext}(x)=\hat{\xi}_{1}(x)-\hat{\xi}_{0}(x)$. In our synthetic
experiments, we note that fitting all three conditional parameters of the
Generalized Extreme Value (GEV) distribution can be numerically unstable.
Given that our primary objective is to investigate the variability of the
shape parameter, we assume that the location and scale parameters are known.
Our focus is to estimate the shape parameter as a function of the covariates,
in order to accurately estimate the CETE and the ETE. It is crucial to
highlight that the shape parameter is our parameter of interest since it
characterizes the heaviness of the tail—a feature that is not encapsulated by
pointwise information. This characteristic cannot be adequately described by
the moments of the distribution or by its location and scale parameters alone.

Figure 2: Convergence of different Extreme Treatment Effect (ETE) estimators
over increasing sample sizes, showcasing the performance of the Naive
Estimator and Proposed Estimator compared to the True ETE. The error bars
represent the variability in the estimates across different simulations.

##  5 Experiments

Figure 3: From left to right, the figures represent the estimation of the
individual tail parameter under the control condition, the tail parameter
under the treatment condition, and the Conditional Extreme Treatment Effect
(CETE) using our proposed method.

In this section, we demonstrate the empirical efficacy of the proposed
framework. First, to verify the estimation procedure for conditional GEV
distributions, we evaluate the proposed algorithm in a synthetic 1D scenario.
Subsequently, we apply our approach to several semi-synthetic high-dimensional
datasets.

Baseline Treatment Effect Models. The proposed framework can be used in
conjunction with existing CATE estimation methods. As such, we consider a few
models used in the literature and adapt them to the proposed framework for
estimating the CETE. The $s$-learner method treats the treatment variable as
another covariate and builds a neural network $f(x,t)$ to estimate the
potential outcomes and then estimates the ITE as
$\hat{\tau}(x)=f(x,1)-f(x,0)$. The $t$-learner builds two neural networks, one
for each treatment group $f_{1}(x)$ and $f_{0}(x)$, and the individual outcome
is given by $\hat{\tau}(x)=f_{1}(x)-f_{0}(x)$ [Künzel et al., 2019]. The
TARNet model is a framework for estimating conditional average treatment
effects with counterfactual balancing [Shalit et al., 2017]. It consists of a
pair of functions, $(\Phi,h)$, where $\Phi$ is a representation learning
function and $h$ is a function that learns the two potential outcomes in the
representation space. CATE is then estimated as
$\hat{\tau}(x)=h(\Phi(x),1)-h(\Phi(x),0)$ The TARNet model minimizes the
objective $\mathcal{L}(\Phi,h)$, which balances the model’s performance on
factual data and the similarity of the representations in the latent space
using the integral probability metric (IPM). The balancing weight, $\alpha$,
controls the trade-off between the similarity and model performance. When
$\alpha=0$ the model is denoted by TARNet, when $\alpha>0$ the model is known
as counterfactual regression (CFR). It corresponds to the 1-Wasserstein
distance for 1-Lipschitz functions Villani [2009].

###  5.1 Synthetic Data Experiments

We commence by validating our proposed method through an experiment involving
synthetic data. We generate $5000$ samples for $Y_{0}|X=x$ from a Generalized
Extreme Value (GEV) distribution, $\text{GEV}(\exp(x),1,\xi_{0}(x))$, where
$\xi_{0}(x)=1+|x|^{0.1}$, and for $Y_{1}|X=x$ from
$\text{GEV}(x^{2},1,\xi_{1}(x))$, with $\xi_{1}(x)=\log(1.1+x^{2})$. The
feature variable $X$ is sampled from a normal distribution,
$\mathcal{N}(0,1)$, and the propensity score is set to $0.3$ for $X>0$ and
$0.7$ otherwise. We evaluate the model on data sampled from $\mathcal{N}(0,1)$
and with a propensity score $p(t=1|x)=\sigma(50x^{2}-5)$, where $\sigma$
denotes the sigmoid function.

The results for this experiment are presented in Figure 3. We observe that the
model accurately learns the shape parameter as a function of the covariates,
which leads to precise estimation of the Conditional Extreme Treatment Effect
(CETE). We then use this trained neural network to estimate the Extreme
Treatment Effect (ETE) by sampling from the learned GEV distribution for
potential outcomes, since we have learned the shape parameter. We fit a GEV to
the marginal samples of the outcomes $\\{\hat{z}^{(j)}_{t}\\}_{j=1}^{N}$. We
also evaluate the naive estimator by fitting a GEV to the observed treatment
outcomes $\\{z^{(j)}_{1}\\}_{j=1}^{N_{1}}$ and
$\\{z^{(j)}_{0}\\}_{j=1}^{N_{0}}$, where $N_{1}$ and $N_{0}$ are the
cardinalities of the treatment and control groups, respectively. We observe
that our estimator exhibits unbiased behavior, unlike the naive estimator. The
convergence behavior of both the naive and the proposed ETE estimators are
presented in Figure 2.

Table 1: Comparison of the CETE performance of our proposed estimator. The error bars are calculated for various random seeds. Model | IHDP (original) | IHDP (Fréchet) | IHDP (Weibull)
---|---|---|---
S-Learner | $0.048\pm 0.028$ | $0.396\pm 0.060$ | $0.307\pm 0.070$
T-Learner | $0.188\pm 0.053$ | $0.634\pm 0.307$ | $0.345\pm 0.004$
TARNet | $0.119\pm 0.071$ | $0.223\pm 0.105$ | $0.278\pm 0.172$
CFR(Wass) | $0.069\pm 0.004$ | $0.211\pm 0.104$ | $0.261\pm 0.061$

###  5.2 Semi-synthetic data experiments

Having validated the proposed method in estimating conditional GEVs, we now
proceed to apply this in the causal inference setting to compute the ETE and
the CETE.

Dataset descriptions. We present a representative family of semi-synthetic
datasets tailored for estimating the CETE. Specifically, we use the Infant
Health and Development Program (IHDP) dataset [Ramey et al., 1992] and two
variants corresponding to different MDAs. Both datasets consider potential
outcomes generated based on known functions of the observed covariates. The
selection bias is the same in the three datasets, provided by the real-world
experiment, and we produce the potential outcomes following different GEVs. We
provide more detailed descriptions of these datasets in Appendix F.

##### CETE results.

We present the performance for CETE estimation in Table 1. The results suggest
that the CFR(Wass) model performs the best for the Fréchet and the Weibull
cases while the SNet performs best for the IHDP original case (the Gumbel
case).

##### ETE results.

We consider the performance of our proposed ETE estimator compared to that of
a naive estimator, which is based on the block maxima method applied directly
to the potential outcomes. The naive estimator is then given by:
$\stackunder[1.2pt]{$\tau$}{\rule{3.44444pt}{0.32289pt}}_{\text{naive}}=\stackunder[1.2pt]{$\xi$}{\rule{3.44444pt}{0.32289pt}}_{1,\text{naive}}-\stackunder[1.2pt]{$\xi$}{\rule{3.44444pt}{0.32289pt}}_{0,\text{naive}}$,
where $\xi_{t,\text{naive}}$ correspond to the estimated shape parameter for
the MDA of the conditional outcome $Y|T=t$. It is important to note that the
conditional outcome $Y|T=t$ is distinct from the potential outcome $Y_{t}$.

Dataset |  $\varepsilon_{\text{ETE}}$(Naive) |  $\varepsilon_{\text{ETE}}$(Ours)
---|---|---
IHDP (original) | $0.094$ | $0.001$
IHDP (Fréchet) | $0.346$ | $0.095$
IHDP (Weibull) | $0.931$ | $0.827$
Table 2: ETE Performance comparison under the naive and proposed estimators
across different datasets using the CFR(Wass) model. Table 3: Overview of the
causal inference datasets for ETE and CETE Estimation.

Name | Type | IPO MDA
---|---|---
IHDP (original) | Semi-Synthetic | Gumbel
IHDP-F | Semi-Synthetic | Fréchet
IHDP-W | Semi-Synthetic | Weibull

IPO MDA: Individual Potential Outcomes Maximum Domain of Attraction.

##  6 Conclusion

In this work, we propose a new framework for estimating extreme treatment
effects. We define new statistical quantities CETE and ETE capturing extreme
treatment effects. We describe conditions to guarantee the identifiability and
the consistency of the proposed estimators. To circumvent the lack of tail
data, we propose a practical approach to sample from the MDAs and introduce a
neural network approach for learning the GEV parameters. Our contributions
offer new insights into the estimation of extreme treatment effects and
provide a basis for further research in this area. Additionally, further
analysis in developing a computational method leveraging the point process
perspective should be considered.

##### Limitations

The main limitation of the method is we require regularity on the space of
conditional outcomes to be able to use our block maxima surrogate.
Circumventing this requires multiple outcomes per individual, which is often
impossible in practice. Understanding any additional structure that can be
used for the block maxima is a direction for future research to address this
limitation.

#### Acknowledgments

This work was supported in part by the National Science Foundation (NSF) under
the National AI Institute for Edge Computing Leveraging Next Generation
Wireless Networks Grant #2112562. Ali Hasan was supported in part by the Air
Force Office of Scientific Research under award number FA9550-20-1-0397.

## References

  * Abadie and Cattaneo [2018] Alberto Abadie and Matias D Cattaneo.  Econometric methods for program evaluation.  _Annual Review of Economics_ , 10:465–503, 2018.
  * Alaa and Van Der Schaar [2017] Ahmed M Alaa and Mihaela Van Der Schaar.  Bayesian inference of individualized treatment effects using multi-task gaussian processes.  _Advances in neural information processing systems_ , 30, 2017.
  * Bingham et al. [1989] Nicholas H Bingham, Charles M Goldie, and Jef L Teugels.  _Regular variation_.  Number 27. Cambridge university press, 1989.
  * Chernozhukov and Hansen [2005] Victor Chernozhukov and Christian Hansen.  An iv model of quantile treatment effects.  _Econometrica_ , 73(1):245–261, 2005.
  * Chernozhukov et al. [2013] Victor Chernozhukov, Iván Fernández-Val, and Blaise Melly.  Inference on counterfactual distributions.  _Econometrica_ , 81(6):2205–2268, 2013.
  * Coles et al. [2001] Stuart Coles, Joanna Bawa, Lesley Trenner, and Pat Dorazio.  _An introduction to statistical modeling of extreme values_ , volume 208.  Springer, 2001.
  * Davison et al. [2012] A. C. Davison, S. A. Padoan, and M. Ribatet.  Statistical Modeling of Spatial Extremes.  _Statistical Science_ , 27(2):161 – 186, 2012\.
  * Díaz [2017] Iván Díaz.  Efficient estimation of quantiles in missing data models.  _Journal of Statistical Planning and Inference_ , 190:39–51, 2017.
  * Firpo [2007] Sergio Firpo.  Efficient semiparametric estimation of quantile treatment effects.  _Econometrica_ , 75(1):259–276, 2007.
  * Fisher and Tippett [1928] Ronald Aylmer Fisher and Leonard Henry Caleb Tippett.  Limiting forms of the frequency distribution of the largest or smallest member of a sample.  In _Mathematical proceedings of the Cambridge philosophical society_ , volume 24, pages 180–190. Cambridge University Press, 1928.
  * Frölich and Melly [2013] Markus Frölich and Blaise Melly.  Unconditional quantile treatment effects under endogeneity.  _Journal of Business & Economic Statistics_, 31(3):346–357, 2013.
  * Galambos and Seneta [1973] Janos Galambos and Eugene Seneta.  Regularly varying sequences.  _Proceedings of the American Mathematical Society_ , 41(1):110–116, 1973.
  * Giessing and Wang [2021] Alexander Giessing and Jingshen Wang.  Debiased inference on heterogeneous quantile treatment effects with regression rank-scores.  _arXiv preprint arXiv:2102.01753_ , 2021.
  * Gnedenko [1943] Boris Gnedenko.  Sur la distribution limite du terme maximum d’une serie aleatoire.  _Annals of mathematics_ , pages 423–453, 1943.
  * Haan and Ferreira [2006] Laurens Haan and Ana Ferreira.  _Extreme value theory: an introduction_ , volume 3.  Springer, 2006.
  * Hill [2011] Jennifer L Hill.  Bayesian nonparametric modeling for causal inference.  _Journal of Computational and Graphical Statistics_ , 20(1):217–240, 2011.
  * Holland [1986] Paul W Holland.  Statistics and causal inference.  _Journal of the American statistical Association_ , 81(396):945–960, 1986.
  * Huang et al. [2022] Wei Huang, Shuo Li, and Liuhua Peng.  Extreme continuous treatment effects: Measures, estimation and inference.  _arXiv preprint arXiv:2209.00246_ , 2022.
  * Imbens [2004] Guido W Imbens.  Nonparametric estimation of average treatment effects under exogeneity: A review.  _Review of Economics and statistics_ , 86(1):4–29, 2004.
  * Imbens and Rubin [2015] Guido W Imbens and Donald B Rubin.  _Causal inference in statistics, social, and biomedical sciences_.  Cambridge University Press, 2015.
  * Johansson et al. [2016] Fredrik Johansson, Uri Shalit, and David Sontag.  Learning representations for counterfactual inference.  In _International conference on machine learning_ , pages 3020–3029. PMLR, 2016.
  * Kallus [2023] Nathan Kallus.  Treatment effect risk: Bounds and inference.  _Management Science_ , 69(8):4579–4590, 2023\.
  * Keele [2015] Luke Keele.  The statistics of causal inference: A view from political methodology.  _Political Analysis_ , 23(3):313–335, 2015.
  * Kennedy [2020] Edward H Kennedy.  Towards optimal doubly robust estimation of heterogeneous causal effects.  _arXiv preprint arXiv:2004.14497_ , 2020.
  * Kennedy et al. [2021] Edward H Kennedy, Sivaraman Balakrishnan, and Larry Wasserman.  Semiparametric counterfactual density estimation.  _arXiv preprint arXiv:2102.12034_ , 2021.
  * Künzel et al. [2019] Sören R Künzel, Jasjeet S Sekhon, Peter J Bickel, and Bin Yu.  Metalearners for estimating heterogeneous treatment effects using machine learning.  _Proceedings of the national academy of sciences_ , 116(10):4156–4165, 2019.
  * Morgan and Winship [2015] Stephen L Morgan and Christopher Winship.  _Counterfactuals and causal inference_.  Cambridge University Press, 2015.
  * Neyman [1923] Jerzy Neyman.  Sur les applications de la théorie des probabilités aux experiences agricoles: Essai des principes.  _Roczniki Nauk Rolniczych_ , 10(1):1–51, 1923\.
  * Pearl [2009] Judea Pearl.  _Causality: Models, Reasoning and Inference_.  Cambridge University Press, USA, 2nd edition, 2009.  ISBN 052189560X.
  * Ramey et al. [1992] Craig T Ramey, Donna M Bryant, Barbara H Wasik, Joseph J Sparling, Kaye H Fendt, and Lisa M La Vange.  Infant health and development program for low birth weight, premature infants: Program elements, family participation, and child intelligence.  _Pediatrics_ , 89(3):454–465, 1992.
  * Resnick [1988] Sidney I Resnick.  Association and multivariate extreme value distributions.  _Australian Journal of Statistics_ , 30(1):261–271, 1988.
  * Resnick [2007] Sidney I Resnick.  _Heavy-tail phenomena: probabilistic and statistical modeling_.  Springer Science & Business Media, 2007.
  * Roncalli [2020] Thierry Roncalli.  _Handbook of financial risk management_.  CRC press, 2020.
  * Rothman and Greenland [2005] Kenneth J Rothman and Sander Greenland.  Causation and causal inference in epidemiology.  _American journal of public health_ , 95(S1):S144–S150, 2005.
  * Rubin [1974] Donald B Rubin.  Estimating causal effects of treatments in randomized and nonrandomized studies.  _Journal of educational Psychology_ , 66(5):688, 1974.
  * Shalit et al. [2017] Uri Shalit, Fredrik D Johansson, and David Sontag.  Estimating individual treatment effect: generalization bounds and algorithms.  In _International Conference on Machine Learning_ , pages 3076–3085. PMLR, 2017.
  * Villani [2009] Cédric Villani.  _Optimal transport: old and new_ , volume 338.  Springer, 2009.
  * Wager and Athey [2018] Stefan Wager and Susan Athey.  Estimation and inference of heterogeneous treatment effects using random forests.  _Journal of the American Statistical Association_ , 113(523):1228–1242, 2018.
  * Zhang [2018] Yichong Zhang.  Extremal quantile treatment effects.  _The Annals of Statistics_ , 46(6B):3707 – 3740, 2018.
  * Zhang et al. [2012] Zhiwei Zhang, Zhen Chen, James F Troendle, and Jun Zhang.  Causal inference on quantiles with an obstetric application.  _Biometrics_ , 68(3):697–706, 2012.

Appendix

##  Appendix A Simulating GEV Potential Outcomes

In this illustration, we model potential outcomes under control and treatment
scenarios for a population of $n=5$ individuals.

Each individual $k\in\\{1,\ldots,n\\}$ is associated with two potential
outcomes. For the control scenario, the potential outcome for individual $k$
follows a Generalized Extreme Value (GEV) distribution with a location
parameter $\mu^{0}_{k}=5+k$, a scale parameter $\sigma^{0}_{k}=1$, and a shape
parameter $\xi^{0}_{k}=0.5$. For the treatment scenario, the potential outcome
for individual $k$ follows a GEV distribution with a location parameter
$\mu^{1}_{k}=6+k$, a scale parameter $\sigma^{1}_{k}=1.5$, and a shape
parameter $\xi^{1}_{k}=1.5+0.1k$.

Figure 4 illustrates the distribution of the potential outcomes for the $n=5$
individuals, each with distinct location and shape parameters. In our work, we
study the differences in the $\xi$ parameter, which, as the figure shows,
dictates how the tail of the distribution behaves.

Figure 4:  Different conditional outcomes for different treatment groups. The
goal of our method is to estimate the difference in the shape parameters of
these limiting GEVs.

##  Appendix B Theoretical Results

In this section, we begin by providing proofs for the remarks presented in the
paper, which support the proposed ETE and CETE definitions. This is followed
by a discussion of other potential definitions for the extreme treatment
effects. Subsequently, we present proofs for the causal inference results
presented in the main text. Finally, we provide a theoretical analysis of the
$\varepsilon$-max sampler algorithm utilized in our empirical settings.

###  B.1 Proof of Regularly Varying Remark

######  Remark 3.3.

Let
$\stackunder[1.2pt]{$\xi$}{\rule{3.44444pt}{0.32289pt}}_{0},\stackunder[1.2pt]{$\xi$}{\rule{3.44444pt}{0.32289pt}}_{1}>0$.
Let $S_{1}(Y_{1}),S_{0}(Y_{0})$ be the survival distributions of the potential
outcome variables with positive support for the treatment and control groups
respectively. Then if
$\stackunder[1.2pt]{$\tau$}{\rule{3.44444pt}{0.32289pt}}_{ext}>0$, then
${S_{0}(Y_{0})}/{S_{1}(Y_{1})}\in
RV(\stackunder[1.2pt]{$\xi$}{\rule{3.44444pt}{0.32289pt}}_{1}-\stackunder[1.2pt]{$\xi$}{\rule{3.44444pt}{0.32289pt}}_{0})$.
Moreover, $S_{1}(Y_{1})<S_{0}(Y_{0})$ for $Y_{t}\gg 1$.

###### Proof.

Define $F_{0},F_{1}$ as the CDFs of the outcome variables for the control and
treatment groups respectively. Then $S_{0}=1-F_{0}$ and $S_{1}=1-F_{1}$.
$F_{0},F_{1}$ are also assumed to be regularly varying with indices of regular
variation as $\xi_{0},\xi_{1}$. By the Karamata characterization theorem, we
can write $F_{0},F_{1}$ in terms of a slowly-varying function, that is, a
function where $\xi=0$. We will call this function $L(Y)$ and then write
$S_{0}=Y^{-\xi_{0}}L_{0}(Y)$ and $S_{1}=Y^{-\xi_{1}}L_{1}(Y)$. As
$Y\to\infty$, the polynomial component dominates, and we can compare the ratio
of the polynomial components. Taking the ratio, we get

$\lim_{Y\to\infty}\frac{S_{0}}{S_{1}}\sim\frac{Y^{-\xi_{0}}}{Y^{-\xi_{1}}}=Y^{\xi_{1}-\xi_{0}}.$

Thus, the difference in the tail indices can be seen as the ratio of the
survival functions for regions deep in the tail. ∎

###  B.2 Proof of Point Process Remark

######  Remark 3.4 (Expected Number of Events).

Consider a causal inference task with potential outcomes variable $Y_{t}$,
$t\in\\{0,1\\}$, that $Y_{t}$ is in the maximum domain of attraction of a GEV
with rate $\xi$, and ETE given by
$\stackunder[1.2pt]{$\tau$}{\rule{3.44444pt}{0.32289pt}}_{ext}>0$. Suppose
that $\mathcal{R}\subset\mathbb{R}$ is the Borel set of interest that should
be quantified for the causal inference task with $|\inf\mathcal{R}|\gg 0$.
Then, $\mathbb{E}[N_{Y_{1}}(\mathcal{R})]>\mathbb{E}[N_{Y_{0}}(\mathcal{R})]$
where
$N_{Y_{t}}(\mathcal{R})=\sum_{i=1}^{\infty}\delta_{Y_{t}^{(i)}\in\mathcal{R}}$
for $t\in\\{0,1\\}$ is a counting process over outcomes on some Borel subset
of $\mathbb{R}$.

###### Proof.

From our assumptions on $\mathcal{R}$, we suppose that the GEV is the
appropriate description of the events. Additionally, suppose that both $Y_{0}$
and $Y_{1}$ have the same limiting sequence $a_{n},b_{n}$ and by EVT

$\mathbb{E}[N_{Y_{0}}(\mathcal{R})]=\left[1+\xi_{0}c_{\mathcal{R}}\right]^{-1/\xi_{0}}$

and

$\mathbb{E}[N_{Y_{1}}(\mathcal{R})]=\left[1+\xi_{1}c_{\mathcal{R}}\right]^{-1/\xi_{1}}.$

Since $\xi_{1}>\xi_{0}$,

$\displaystyle 1+\xi_{1}c_{\mathcal{R}}$
$\displaystyle>1+\xi_{0}c_{\mathcal{R}}$
$\displaystyle(1+\xi_{1}c_{\mathcal{R}})^{1/\xi_{1}}$
$\displaystyle>(1+\xi_{0}c_{\mathcal{R}})^{1/\xi_{1}}$
$\displaystyle(1+\xi_{1}c_{\mathcal{R}})^{1/\xi_{1}}$
$\displaystyle<(1+\xi_{0}c_{\mathcal{R}})^{1/\xi_{0}}$
$\displaystyle(1+\xi_{1}c_{\mathcal{R}})^{-1/\xi_{1}}$
$\displaystyle>(1+\xi_{0}c_{\mathcal{R}})^{-1/\xi_{0}}$

with $c_{\mathcal{R}}\in\mathbb{R}_{+}\cup\\{0\\}$. ∎

###  B.3 Possible Definitions for Extreme Treatment Effects

In this section, we note that while we motivate our definition for the extreme
treatment effect from an intuitive perspective, a point process perspective,
and from a regular variation perspective, it is possible to define other
functions of the shape parameter $\Psi(\xi_{1}(x),\xi_{0}(x))$ that capture
other properties. Our presented methods can be extended to these general
functions, as the main problem in causal inference lies in estimating
$\xi_{1}(x)$ and $\xi_{0}(x)$. Moreover, the choice of the difference between
the shape parameters guarantees that the statistical properties of the
estimators of shape parameters (e.g., expectation and variance) can be easily
extended to the estimator of the extreme treatment effects, which may not be
the case for other functions such as the ratio of the shape parameters. It is
important to highlight that, at times, quantifying the proportional change of
the shape parameter can be essential (thus, taking the ratio would be more
justified). In contrast, in some situations, the magnitude of change itself is
significant regardless of the initial value (for instance, if having a heavier
tail above $2$ is critical and the control shape parameter is $0$, taking the
difference provides more information). Therefore, different definitions might
be more relevant depending on the context

###  B.4 Causal Inference Proofs

We restate the main statements and then state the proofs.

######  Proposition 3.6.

If unconfoundedness holds then tail unconfoundedness holds.

###### Proof.

Assuming that unconfoundedness holds, and letting $t\in\mathcal{T}$, we have

$Y_{t}\perp\\!\\!\\!\\!\perp T\mid X$

Therefore,

$P(Y_{t}\mid X,T)=P(Y_{t}\mid X).$

For $\\{Y^{(k)}\\}_{k=1}^{m}$ i.i.d copies of $Y$, we have that

$\forall k\in\\{1,\ldots,m\\},\;P(Y^{(k)}_{t}\mid X,T)=P(Y^{(k)}_{t}\mid X)$

By taking the limit on both sides, we have

$\lim_{m\rightarrow\infty}{P(Z_{t}^{m}\mid
X,T)=\lim_{m\rightarrow\infty}P(Z_{t}^{m}\mid X)}$

which implies tail unconfoundedness holds. ∎

######  Proposition 3.9 (Identifiability).

Under Assumptions 2.2, 2.3, and 3.5, the CETE and ETE estimators defined in
Definition 3.7 and Definition 3.8 respectively are causally identifiable.
Moreover, for $t\in\\{0,1\\}$, the proposed estimator satisfies

$\hat{\xi}_{t}=\arg\max_{\theta\in\Theta}\sum_{j=1}^{N}\log p\left(Z^{(j)}\mid
g_{\theta}(x^{(j)}),t^{(j)}\right)\mathbbm{1}_{t^{(j)}=t}.$

###### Proof.

By definition, we have

$\displaystyle\hat{\theta}_{t}=$
$\displaystyle\arg\max_{\theta}\lim_{m\to\infty}\sum_{j=1}^{N}\log
P\left(Z^{m}_{t}(x^{(j)})\;\big{|}\;\theta\left(x^{(j)}\right)\right)$
$\displaystyle=$
$\displaystyle\arg\max_{\theta}\sum_{j=1}^{N}\lim_{m\to\infty}\log
P\left(Z^{m}_{t}(x^{(j)})\;\big{|}\;\theta\left(x^{(j)}\right)\right)$

By tail unconfoundedness, we have

$\lim_{m\to\infty}P\left(Z^{m}_{t}(x^{(j)})\;\big{|}\;\theta\left(x^{(j)}\right)\right)=\lim_{m\to\infty}P\left(Z^{m}_{t}(x^{(j)})\;\big{|}\;\theta\left(x^{(j)}\right),T=t\right)$

From the consistency assumption, it follows that

$\lim_{m\to\infty}P\left(Z^{m}_{t}(x^{(j)})\;\big{|}\;\theta\left(x^{(j)}\right),T=t\right)=\lim_{m\to\infty}P\left(Z^{m}(x^{(j)})\;\big{|}\;\theta\left(x^{(j)}\right),T=t\right)$

Notice that the right-hand-side quantity does not involve any potential
outcomes. We then have

$\lim_{m\to\infty}P\left(Z^{m}_{t}(x^{(j)})\;\big{|}\;\theta\left(x^{(j)}\right)\right)=\lim_{m\to\infty}P\left(Z^{m}(x^{(j)})\;\big{|}\;\theta\left(x^{(j)}\right),T=t\right)$

Hence, the estimator is identifiable. ∎

######  Remark 3.10.

The ETE estimator defined in Definition 3.8 is consistent and asymptotically
normal.

###### Proof.

Given that the MLE estimator for the shape parameters is consistent, we can
apply Slutsky’s theorem to deduce that the CETE estimator is also consistent,
for a fixed individual as his potential outcomes observations go to infinity.
Additionally, because the MLE estimator for the shape parameters is
asymptotically normal, we can use the Delta method to conclude that the ETE
estimator is asymptotically normal. ∎

##  Appendix C Theoretical Analysis for $\varepsilon$-max-sampler

The main idea of this approach is to consider a ball around each individual
outcome and then use the maxima of all outcomes within this ball as an
estimate of data in the tails. We prove that this converges to the true maxima
under suitable regularity conditions. For treatment groups $t\in\\{0,1\\}$,
the conditional potential outcome variable
$\\{(Y|X=x)=Y_{t}(x),\;x\in\mathcal{X}\\}$ forms a stochastic process over
$\mathcal{X}$. For every $x\in\mathcal{X}$, we assume that
$\left(Y^{(1)}_{t}(x),Y^{(2)}_{t}(x),\ldots,Y^{(m)}_{i}(x)\right)$ are iid
random variables with the same distribution as $Y(x)$. We define the max-
sampler: $Z^{m}_{t}(x)=\max_{1\leq k\leq m}Y^{(k)}_{t}(x),\quad
x\in\mathcal{X}.$ Let $\varepsilon>0$, and let
$B_{\varepsilon}(x)=\\{x^{\prime}\mid
x^{\prime}\in\mathbb{R}^{d}\;\text{and}\;\|x-x^{\prime}\|<\varepsilon\\}$ be
the open ball around $x$ with radius $\varepsilon$. Let
$B^{m}_{\varepsilon}(x)$ be a subset of $m$ points of $B_{\varepsilon}(x)$. We
define the $\varepsilon$-max-sampler as

$\hat{Z}^{m}_{t,\epsilon}(x)=\max_{x^{\prime}\in
B^{m}_{\epsilon}(x)}Y_{t}(x^{\prime})$ (2)

We next study the theoretical relationships between $Z^{m}_{t}(x)$ and
$\hat{Z}^{m}_{t,\varepsilon}(x)$ and determine the conditions under which
$\hat{Z}^{m}_{t,\varepsilon}(x)$ converges to $Z^{m}_{t}(x)$.

######  Proposition C.1 (Convergence of $\varepsilon$-max-sampler).

Let $t\in\\{0,1\\}$ and $m\in\mathbb{N}$. Assume that for every
$x,x^{\prime}\in\mathcal{X}$, such that $x\neq x^{\prime}$, we have that
$Y(x)\perp\\!\\!\\!\\!\perp Y(x^{\prime})$ and $x\mapsto P(Y(x)\leq y)$ is a
continuous function of $x$. Then for every limit point, $x\in\mathcal{X}$, the
following convergence holds in distribution

$\forall
x\in\mathcal{X},\;\hat{Z}^{m}_{t,\varepsilon}(x)\underset{\varepsilon\to
0}{\longrightarrow}Z^{m}_{t}(x)$

###### Proof.

Denote by $F(\cdot\mid x)$ the CDF of $\hat{Z}^{m}_{t,\varepsilon}(x)$ and let
$x_{1},\ldots,x_{n}$ be $n$ elements of the $\varepsilon$-ball. We have that

$\displaystyle F(z|x)$ $\displaystyle=P(\hat{Z}^{m}_{t,\varepsilon}(x)\leq z)$
$\displaystyle=P(Y(x_{1})\leq z,\ldots,Y(x_{n})\leq z)$
$\displaystyle=\Pi_{k=1}^{n}P(Y(x_{k})\leq z)$

we have that when $\epsilon\to 0$ all the points in the $\varepsilon-$ball
$x_{k}\to x$, therefore,

$\Pi_{k=1}^{n}P(Y(x_{k})\leq z)\to\Pi_{k=1}^{n}P(Y^{(k)}(x)\leq
z),\;\text{as}\;\varepsilon\to 0.$

We then have that

$\Pi_{k=1}^{n}P(Y^{(k)}(x)\leq z)=P(Y^{(1)}(x)\leq z,\ldots,Y^{(n)}(x)\leq z)$

and the right-hand side corresponds to the CDF of $Z^{m}_{t}(x)$. ∎

We note that the first independence assumption is often referred to as the no
interference assumption between the units (see also the stable unit treatment
value assumption [Imbens and Rubin, 2015, Keele, 2015]). From this, we can
bound the difference between the $\varepsilon$-max-sampler estimate and the
true potential outcome under a Lipschitz assumption on the potential outcomes:

######  Proposition C.2.

Let $t\in\\{0,1\\}$. Assume that the potential outcomes are positive and that
$x\rightarrow Y_{t}(x)$ is K-Lipschitz almost surely, then

$\forall
n\in\mathbb{N},\forall\varepsilon>0,\;\left\|\hat{Z}^{m}_{t,\varepsilon}(x)-Z^{m}_{t}(x)\right\|\leq
K\varepsilon\;\;\text{a.s.}$

The theorem states that the error grows at most at the rate of the Lipschitz
constant.

###### Proof.

Let $n\in\mathbb{N}$ and let $\varepsilon>0$. We have that

$\forall x,x^{\prime}\in\mathcal{X},\;\left\|Y(x)-Y(x^{\prime})\right\|\leq
K\left\|x-x^{\prime}\right\|$

For the $\varepsilon-$ball around $x$ we have that

$\forall x^{\prime}\in
B_{\varepsilon}(x),\;\left\|Y(x)-Y(x^{\prime})\right\|\leq
K\left\|x-x^{\prime}\right\|.$

By taking the max over the right-hand side, we have that

$\forall x^{\prime}\in
B_{\varepsilon}(x),\;\left\|Y(x)-Y(x^{\prime})\right\|\leq K\varepsilon.$

Since the potential outcomes are positive we have for $x_{1},\ldots,x_{m}\in
B_{\varepsilon}(x)$

$\left\|\max_{1\leq k\leq m}Y(x_{k})-\max_{1\leq k\leq
m}Y^{k}(x)\right\|\leq\max_{1\leq k\leq m}\left\|Y(x_{k})-Y^{k}(x)\right\|.$

Therefore,

$\left\|\hat{Z}^{m}_{t,\varepsilon}(x)-Z^{m}_{t}(x)\right\|\leq K\varepsilon.$

∎

##  Appendix D Tail Unconfoundednes Examples

In this section, we present two quantitative examples where tail
unconfoundedness is satisfied but unconfoundedness does not hold.

##### Example 1:

Let the mass function $P(Y_{1}\mid T=t,X)$ be a Bernoulli distribution with a
non-binary treatment parameter $t\in(0,1)$. To compute the probability that
the maximum is equal to 1, we write

$\displaystyle P(Z_{1}^{n}=1\mid T=t,X)$
$\displaystyle=1-P(Y^{1}_{1}=0,\ldots,Y^{n}_{1}=0\mid T=t,X)$
$\displaystyle=1-P^{n}(Y_{1}\mid T=t,X)$ $\displaystyle=1-t^{n}.$

This is independent of $t$ as $n\to\infty$ (since $t\in(0,1)$).

##### Example 2:

Consider the density function $p(y_{1}\mid x,t)$ that follows a Beta
distribution $\mathcal{B}(\alpha,\beta)$ with parameters $\alpha=x$ and
$\beta=t+1$ for $x,t>0$. The domain of attraction of this distribution is
$\text{GEV}(0,1,\xi(x))$ with $\xi(x)=-x^{-1}$ with the following normalizing
constants Roncalli [2020]:

$\displaystyle b_{m}(x)=1\quad\text{and}\quad
a_{m}(x)=\left(\frac{m\Gamma(x+t+1)}{\Gamma(x)\Gamma(t+1)}\right)^{-1/(t+1)}.$

We see that as $m\to\infty$, $a_{m}\to 0$. Then, the limiting distribution of
$Z_{1}^{n}$ only depends on $x$ and does not depend on $t$, implying that the
tail unconfoundedness holds.

##  Appendix E Algorithm

We describe the algorithm for the CETE estimation in algorithm 1.

1

Data: $S=\\{\left(x^{(i)},t^{(i)},y^{(i)}\right)\\}_{i=1}^{n}$

Input: A neural network
$N_{\theta}=\\{(\hat{\mu}_{t},\hat{\sigma}_{t},\hat{\xi}_{t})\\}_{t=0}^{1}$ :
$s$-Learner, $t$-Learner, TARNet, or CFR

Output: A trained neural network $N_{\theta^{*}}$

2

3 Function  _MaxSampler(_$S,m$_)_:

4 Initialize an empty set $S^{m}$

5 Number of clusters:

6 $K=n\;\text{div}\;m$

7 Run $K$-means to and cluster the data into $K$ clusters.

8 for  _$i=1,2,\ldots,K$_ do

9 Take $S_{i}=\\{j\mid(x^{(j)},t^{(j)},y^{(j)})\;\text{in
the}\;i^{th}\;\text{cluster}\\}$

10

11 Select $s=\arg\underset{j\in S_{i}}{\text{max}}y^{(j)}$

12

13 Add $(x^{(s)},t^{(s)},y^{(s)})$ to $S^{m}$

14

15 return $S^{m}$

16

17 Function  _Loss(_$S^{m},N_{\theta}$_)_:

18 For each $(x^{(i)},t^{(i)},y^{(i)})\in S^{m}$ the negative conditional log-
likelihood is

19 $\begin{aligned}
L_{i}=&-\log(P(y^{(i)}\mid\mu_{0}(x^{(i)}),\sigma_{0}(x^{(i)}),\xi_{0}(x^{(i)})))\mathbf{1}_{t^{(i)}=0}\\\
&-\log(P(y^{(i)}\mid\mu_{1}(x^{(i)}),\sigma_{1}(x^{(i)}),\xi_{1}(x^{(i)})))\mathbf{1}_{t^{(i)}=1}\end{aligned}$

20 return $\frac{1}{n}\sum_{i=1}^{n}L_{i}$

21 Function  _Main_ :

22 Choose $m$

23 Run the max-sampler

24 $S^{m}=\textnormal{{MaxSampler(}}\textnormal{\emph{S,m}}\textnormal{{)}}$

25 Train the neural network $N_{\theta}$ by minimizing the loss:

26 L = Loss(_$S^{m},N_{\theta}$_)

27

28 return The trained neural network $N_{\theta*}$

29

Algorithm 1 Learning the Parameters of the Domain of Attraction of Individual
Potential Outcomes

##  Appendix F Datasets Descriptions

In this section, we present a detailed description of the datasets used in the
numerical study. A summary is provided in Table 4.

##### IHDP

The IHDP dataset was first introduced by Hill [2011] and it is based on real
covariates available from the Infant Health and Development Program (IHDP),
which studied the effect of development programs on children. The features in
this dataset come from a randomized controlled trial. The potential outcomes
are simulated using Setting B. The dataset consists of $747$ individuals
(specifically, $139$ in the treatment group and $608$ in the control group),
each with $25$ features. The potential outcomes are generated according to the
following distributions for the control group:

$Y_{0}=\exp(\beta^{T}\cdot(X+W))+\eta_{0},\;\;\text{with}\;\eta_{0}\sim\mathcal{N}(0,1)$

and for the treatment group:

$Y_{1}=\beta^{T}(X+W)-\omega+\eta_{1},\;\;\text{with}\;\eta_{1}\sim\mathcal{N}(0,1)$

where $W$ has the same dimension as $X$ with all entries equal $0.5$ and
$\omega=4$. The regression $25$ dimensional coefficient vector $\beta$ is
randomly sampled from a categorical distribution with the support
$(0,0.1,0.2,0.3,0.4)$ and the respective probabilities
$\mu=(0.6,0.1,0.1,0.1,0.1)$. We refer to this dataset as IHDP (Original),
which corresponds to a Gumbel domain of attraction. We also generate two new
versions of IHDP by changing the noise variables of the potential outcomes.
The first version is given by:

$Y_{0}=\exp(\beta^{T}\cdot(X+W))+\eta_{0},\;\;\text{with}\;\eta_{0}\sim\text{GEV}(0,1,0.1+\|X\|_{2}^{1/10})$

and

$Y_{1}=\beta^{T}(X+W)-\omega+\eta_{1},\;\;\text{with}\;\eta_{1}\sim\text{GEV}\left(0,1,\log(1.1+\|X\|_{2})\right)$

We refer to this version as IHDP (Fréchet). The second version is generated
with the following potential outcomes:

$Y_{0}=\exp(\beta^{T}\cdot(X+W))+\eta_{0},\;\;\text{with}\;\eta_{0}\sim\text{GEV}(0,1,-0.1-\|X\|_{2}^{1/10})$

and

$Y_{1}=\beta^{T}(X+W)-\omega+\eta_{1},\;\;\text{with}\;\eta_{1}\sim\text{GEV}\left(0,1,-\log(1.1+\|X\|_{2})\right)$

We refer to this version as IHDP (Weibull).

Table 4: Overview of the causal inference datasets for CETE and ETE
Estimation.

Name | Type | IPO MDA
---|---|---
IHDP (original) | Semi-Synthetic | Gumbel
IHDP (Fréchet) | Semi-Synthetic | Fréchet
IHDP (Weibull) | Semi-Synthetic | Weibull

Individual Potential Outcomes Maximum Domain of Attraction (IPO MDA) if
consistent across all individuals.

##  Appendix G GEV Likelihoods

For reference, we provide the correspondences between the bulk distributions
and tail distributions for our simulated experiments.

####  G.0.1 Gumbel Distribution

The density function of a Gumbel is

$f(y;\mu(x),\beta(x))=\frac{1}{\beta(x)}\exp\left({-\left(\frac{y-\mu(x)}{\beta(x)}+e^{-\frac{y-\mu(x)}{\beta(x)}}\right)}\right)$

with $\mu(x)$ a real function for the location parameter and $\beta(x)$ a
positive real function of the scale.

Then, the conditional negative log-likelihood is

$L(\mu,\beta)=\sum_{i=1}^{n}\left(\frac{y_{i}-\mu(x_{i})}{\beta(x_{i})}+e^{\frac{y_{i}-\mu(x_{i})}{\beta(x_{i})}}\right)+n\log{\beta(x_{i})}.$

####  G.0.2 Fréchet Distribution

The density function of a Fréchet distribution is

$\displaystyle
f(y;\alpha(x),s(x),m(x))=\frac{\alpha(x)}{s(x)}\left(\frac{y-m(x)}{s(x)}\right)^{-1-\alpha(x)}e^{-\left(\frac{y-m(x)}{s(x)}\right)^{-\alpha(x)}}.$

Then, the conditional negative log-likelihood is

$\displaystyle
L(\alpha,s,m)=\sum_{i=1}^{n}\Bigg{[}(1+\alpha(x_{i}))\log\left(\frac{y_{i}-m(x_{i})}{s(x_{i})}\right)+\left(\frac{y_{i}-m(x_{i})}{s(x_{i})}\right)^{-\alpha(x_{i})}+\log\left(\frac{s(x_{i})}{\alpha(x_{i}}\right)\Bigg{]}.$

####  G.0.3 Weibull Distribution

The density function of a Weibull is

$\displaystyle
f(y;\lambda(x),k(x))=\frac{k(x)}{\lambda(x)}\left(\frac{y}{\lambda(x)}\right)^{k(x)-1}\exp{\left(-(y/\lambda(x))^{k(x)}\right)}\mathbbm{1}_{y\geq
0}$

with $k(x)$ and $\lambda(x)$ two positive functions of the shape and the
scale, respectively. The conditional negative log-likelihood is

$\displaystyle
L(\lambda,k)=\sum_{i=1}^{n}\left(\frac{y_{i}}{\lambda(x_{i})}\right)^{k(x_{i})}+k(x_{i})\log(\lambda(x_{i}))-\log(k(x_{i}))+(1-k(x_{i}))\log(y_{i}).$

##  Appendix H Convergence to GEV Distributions

####  H.0.1 Gaussian to Gumbel

Let $Y(x)\sim\mathcal{N}(\mu(x),\sigma(x))$. The limiting GEV is given by

$\frac{Z_{m}(x)-b_{m}(x)}{a_{m}(x)}\underset{m\to+\infty}{\longrightarrow}\text{Gumbel}(0,1)$

where

$b_{m}=\mu(x)+\sigma(x)\left(\sqrt{2\log(m)}-\frac{\log(\log(m))+\log(4\pi)}{2\sqrt{\log(m)}}\right)$

and

$a_{m}=\frac{\sigma(x)}{\sqrt{2\log(m)}}.$

This implies that for large enough $m$,

$Z_{m}(x)\approx c_{m}(x)g+d_{m}(x)$

with $g\sim\text{Gumbel}(0,1)$.

####  H.0.2 Log-Gamma to Fréchet

Let $Y(x)\sim\mathcal{LG}(\alpha(x),\beta(x))$. The limiting GEV is given by

$\frac{Z_{m}(x)-b_{m}(x)}{a_{m}(x)}\underset{m\to+\infty}{\longrightarrow}\text{GEV}(0,1,\beta(x)^{-1})$

where

$a_{m}(x)=0\;\text{and}\;b_{m}(x)=\frac{\left(m\log(m)^{\alpha(x)-1}\right)^{1/\beta(x)}}{\Gamma(\alpha(x))}.$

Then,

$Z_{m}(x)\approx a_{m}(x)\;\text{GEV}(0,1,\xi(x))+b_{m}(x)$

with

$\xi(x)=\beta(x)^{-1}.$

####  H.0.3 Beta to Weibull

Let $Y(x)\sim\mathcal{B}(\alpha(x),\beta(x))$. We have that

$\frac{Z_{m}(x)-b_{m}(x)}{a_{m}(x)}\underset{m\to+\infty}{\longrightarrow}\text{GEV}(0,1,-\alpha(x)^{-1})$

where

$b_{m}(x)=1\;\;\text{and}\;\;a_{m}(x)=\left(\frac{m\Gamma(\alpha(x)+\beta(x))}{\Gamma(\alpha(x))\Gamma(\beta(x)+1)}\right)^{-1/\beta(x)}.$

Then,

$Z_{m}(x)\approx a_{m}(x)\;\text{GEV}(0,1,\xi(x))+b_{m}(x)$

with

$\xi(x)=-\alpha(x)^{-1}.$

(a) Gumbel, $D=1$

(b) Weibull, $D=1$

(c) Frechet, $D=1$

(d) Gumbel, $D=5$

(e) Weibull, $D=5$

(f) Frechet, $D=5$

(g) Gumbel, $D=10$

(h) Weibull, $D=10$

(i) Frechet, $D=10$

Figure 5: Experiments For Conditional GEV Performance as a function of $N$ and
$D$

(a) $D=1$

(b) $D=5$

(c) $D=10$

Figure 6: Experiments for the Gaussian to Gumbel case as a function of K

##  Appendix I Additional Experiments

In this section, we provide additional experiments to validate the proposed
method. We start with verifying our methods for estimating conditional GEVs
and the estimated scale parameter of MDA for the Log-Gamma distribution. Then,
an ablation is performed on the number of the data size $N$ and the
dimensionality $D$. Finally, we present an ablation on the number of clusters
$K$.

####  I.0.1 Estimation from GEV observations

(a) Gumbel distribution with the true location and scale functions
$\mu(x)=x^{2}$ and $\beta(x)=0.25+\frac{|x|}{2}$.

(b) Weibull distribution with the true scale and shape functions
$\lambda(x)=1+x^{2}$ and $\beta(x)=1+|x|$.

(c) Fréchet distribution with the true location and scale functions
$m(x)=e^{x}$, $s(x)=1+|x|$, and $\alpha(x)=1+x^{2}$.

Figure 7: Estimation of distribution parameters from GEV observations for
Gumbel, Weibull, and Fréchet distributions.

We estimate the functional parameters of the GEV given a dataset
$S=\\{(x^{(i)},y^{(i)})\\}_{i=1}^{n}$ where $P(Y\mid X=x)$ follows a GEV
distribution. We consider the cases of the Weibull, Gumbel, and Fréchet
distributions by learning the parameters of each of the distributions. The
full likelihood functions for each case were detailed in Appendix G. A
comparison of the estimated and true functions is given in Figure 7.

##### Conditional Gumbel

We assume that the data are generated according to
$Y(x)\sim\text{Gumbel}(\mu(x),\beta(x))$ with $\mu(x)=x^{2}$ and
$\beta(x)=0.25+\frac{|x|}{2}$. Figure 7(a) shows the performance of our model
on this dataset.

(a) Estimation of parameters of the limiting Gumbel distribution generated by
a Gaussian with $\mu(x)=|x|$ and $\sigma(x)=x^{2}$.

(b) Estimation of parameters of limiting Weibull generated by a beta
distribution with the parameters given by $\alpha(x)=|x|+1$ and
$\beta(x)=x^{2}+1$.

(c) Estimation of parameters of limiting Fréchet distribution generated by a
log-gamma distribution with parameters given by $\alpha(x)=|x|+1$ and
$\beta(x)=x^{2}+1$.

Figure 8: Estimation of limiting GEV distribution parameters from observations
of Gaussian, beta, and log-gamma distributions. All parameters are estimated
by taking the $\max$ of the $K$-means partition of the feature space.

##### Conditional Weibull

We assume that the data are generated according to
$Y(x)\sim\text{Weibull}(\lambda(x),k(x))$ with $\lambda(x)=1+x^{2}$ and
$\beta(x)=1+|x|$. Figure 7(b) illustrates the results of fitting the estimator
to the conditional data.

##### Conditional Fréchet

We assume that the data are generated according to:
$Y(x)\sim\text{Fréchet}(\alpha(x),s(x),m(x))$ where $\alpha(x)=1+x^{2}$,
$s(x)=1+|x|$, and $m(x)=\exp(x)$. Figure 7(c) illustrates the results for
fitting in the Fréchet case.

These experiments validate the estimation procedure in the main text when
observing the extreme data. Next, we consider the case where we observe
samples of the data and need to use the $\varepsilon$-max-sampler to estimate
the extreme observations. We do this by applying the $\varepsilon$-max-sampler
to data from a log-normal distribution and estimate the corresponding Frechét
distribution. The results are presented in Figure 9 for estimating the scale
$\sigma$ parameter.

We also estimated the scale parameter for the Log-Gamma distribution in Figure
9.

Figure 9: Estimation of the limiting scale parameter $\sigma$ of the maximum
of the log-gamma distribution.

####  I.0.2 Estimating conditional limiting GEV distributions

For these experiments, we consider the $\varepsilon$-max-sampler with a total
of $2000$ random samples of the limiting GEV. In Appendix I, we include
detailed experiments on fitting conditional GEV distributions when they are
directly observed. We generate data based on a distribution with a known MDA,
and then compare the recovered GEV distribution parameters to sequences known
for their convergence to the true parameters. We provide full details of the
convergence of these distributions to their corresponding GEV distributions in
Appendix H. The results are presented in Figure 8 where we compare the true
functions to the estimated ones for the different datasets. We additionally
provide an ablation study on the impact of varying the block sizes and the
dimensionality of the data in Appendix I.

Gaussian distribution: Gumbel MDA. We generate $X\sim\mathcal{N}(0,1)$, and
$Y(x)\sim\mathcal{N}(\mu(x),\sigma(x))$, where $\mu(x)=|x|$ and
$\sigma(x)=x^{2}$. The corresponding GEV is the Gumbel distribution. Moreover,
we estimate the limiting conditional GEV parameters with
$\hat{\mu}(x),\hat{\sigma}(x),\hat{\xi}(x)$. Figure 8(a) illustrates this
result where we compare a sequence converging to the ground truth to the
estimated parameters.

Beta distribution: Weibull MDA. We generate samples $X\sim\mathcal{N}(0,1)$
and $Y(x)\sim\mathcal{B}(\alpha(x),\beta(x))$ with $\alpha(x)=|x|+1$ and
$\beta(x)=x^{2}+1$. The corresponding GEV is a Weibull distribution. We
illustrate the estimates of these parameters in Figure 8(b).

Log-Gamma distribution: Fréchet MDA. We generate samples
$X\sim\mathcal{N}(0,1)$ and $Y(x)\sim\mathcal{LG}(\alpha(x),\beta(x))$ with
$\alpha(x)=|x|+1$ and $\beta(x)=x^{2}+1$ where $\mathcal{LG}$ is the log-gamma
distribution. The corresponding GEV is a Fréchet distribution. Figure 8(c)
illustrates the performance of the proposed method. The results for the scale
parameter are presented in Appendix I

##### Ablation on $N$ and $D$

We study how the performance of our method changes when varying the data size
$N$ and the dimensionality $D$, as illustrated in Figure 5. The ablation
suggests that the Weibull and Frechet cases scale when observing large $N$ but
the Gumbel case is more difficult to estimate.

##### Ablation on number of clusters

We study the performance of our method when varying the number of clusters $K$
and the dimensionality $D$ for a fixed number of data points $N=10^{4}$. The
results are illustrated in Figure 6.