Joint Exceedance Probability for Multivariate Distribution I was reading this paper:
https://nhess.copernicus.org/preprints/nhess-2020-28/nhess-2020-28.pdf
Where the notion of "joint exceedance probability" is discussed for the bivariate case. That is:
(1 >  ∩ 2 > )
where  is some threshold (e.g. 90th percentile). This can be estimated by looking at the number of observations for 1 and 2 that exceed both th percentiles relative to the total number of observations. So, say we have 100,000 total observations for the joint CDF of 1 and 2. We define the threshold  as the 90th percentile and count the number of observations in the region where both 1 and 2 > 0.90 as, say 15. We can then (naively) calculate the joint exceedance probability as:
(1 > 0.9 ∩ 2 > 0.9) = 15/100000 = 0.00015 or 0.015%.
My question is: does the bivariate case generalize to higher dimensions? If we have three or more random variables, can we simply say that the joint exceedance probability for the d-dimensional case is:
(1 >  ∩ 2 >  ∩ 3 >  ∩ ... ∩ Xd > )
And use the same idea described above? That is, look for observations in the region where all variables exceed some threshold, count those observations, and divide by the total number of observations?
The answer seems obvious to me (yes)... but maybe I'm missing something.

Yes, it does generalise
In general, for any event $\mathbf{X} \in \mathscr{X}$ involving a random vector $\mathbf{X}$, we can estimate the probability of that event by its empirical estimator based on IID observations.  That is, if $\mathbf{X}_1, \mathbf{X}_2, ..., \mathbf{X}_n$ are IID observations then we can estimate the probability $\mathbb{P}(\mathbf{X} \in \mathscr{X})$ using the empirical estimator:
$$\widehat{\mathbb{P}(\mathbf{X} \in \mathscr{X})} = \frac{1}{n} \sum_{i=1}^n \mathbb{I}(\mathbf{X}_i \in \mathscr{X}).$$
This estimator is unbiased and consistent in the case where the underlying data is IID.
