Multivariate Extreme Value Theory Based Channel Modeling for Ultra-Reliable Communications

Attaining ultra-reliable communication (URC) in fifth-generation (5G) and beyond networks requires deriving statistics of channel in ultra-reliable region by modeling the extreme events. Extreme value theory (EVT) has been previously adopted in channel modeling to characterize the lower tail of received powers in URC systems. In this paper, we propose a multivariate EVT (MEVT)-based channel modeling methodology for tail of the joint distribution of multi-channel by characterizing the multivariate extremes of multiple-input multiple-output (MIMO) system. The proposed approach derives lower tail statistics of received power of each channel by using the generalized Pareto distribution (GPD). Then, tail of the joint distribution is modeled as a function of estimated GPD parameters based on two approaches: logistic distribution, which utilizes logistic distribution to determine dependency factors among the Fréchet transformed tail sequence and obtain a bi-variate extreme value model, and Poisson point process, which estimates probability measure function of the Pickands angular component to model bi-variate extreme values. Finally, validity of the proposed models is assessed by incorporating the mean constraint on probability measure function of Pichanks coordinates. Based on the data collected within the engine compartment of Fiat Linea, we demonstrate the superiority of proposed methodology compared to the conventional extrapolation-based methods in providing the best fit to the multivariate extremes.


I. INTRODUCTION
Ultra-reliability is one of the most important features of the fifth generation (5G) and beyond networks with the goal of providing a reliable and resilient communication infrastructure that can support a wide range of applications and use cases, including industrial automation, self-driving cars, and remote medical procedures [1]- [4].In the current terminology of 5G standards, ultra-reliability is commonly associated with low latency and creates an ultra-reliable low latency communication (URLLC).However, this tight coupling should be relaxed in certain scenarios such as health monitoring or disaster recovery applications where ultra-reliability is essential, but the allowed latency can be larger than the typical 1 milliseconds (ms) [5].Fulfillment of the reliability requirements of 10 −9 -10 −5 packet error rate (PER) in ultra-reliable communication (URC) has been widely used in previous studies [4]- [10].URC necessitates fundamental breakthrough in the derivation of the channel lower tail distribution by using novel statistical tools, such as extreme value theory (EVT) [11], [12], power law expression based on the extrapolation of the channel data [4], [5], and data-driven learning framework [6].Different forms of diversity techniques can be applied in URC to reduce the required signal-to-noise-ratio (SNR) for achieving a certain reliability.Among the diversity techniques, spatial diversity is the most viable solution to achieve high reliability over fading channels by utilizing multiple transmission antennas [13].Consequently, towards achieving ultra-reliability in spatial diversity, the key is to obtain a proper multi-channel modeling approach that characterizes the statistics of multivariate extreme events for the systems with multiple input multiple outputs (MIMO) connections [14].
Ultra-reliability can be achieved through different forms of diversity techniques, including time diversity [15], frequency diversity [16], interface diversity [17], or spatial diversity [4], [18], [19]- [21].In [15], a combination of the stochastic geometry and the queueing analysis is proposed to derive the URLLC achievable ratio through the concept of effective bandwidth, which transforms the constraints of delay and reliability into the constraint of the achievable transmission rate.The authors in [16] explore the diversity aspects of random access schemes where users transmit over multiple orthogonal Rayleigh fading subchannels and are treated as diversity branches by a maximum ratio-combining (MRC) receiver.The authors in [17] consider a trade-off between transmission latency and reliability by allocating coded fragments of the encoded payload message to different interfaces according to their bit rate, latency, and reliability properties.[19] proposes a statistical model by incorporating the time of arrival (TOA), angle of arrival (AOA), and angle of departure (AOD) of each multi-path component to estimate the joint magnitude and phase PDFs of measured data for realistic model parameters.A survey presented in [20] investigates the main properties of massive MIMO channels based on different configurations of an antenna array, which directly affect the channel models and system performance.On the other hand, [21] provides a survey of the most critical concepts in propagation channel modeling for MIMO systems by classifying the channels into physical models focusing on double-directional propagation and analytical models considering the channel impulse response and antenna properties.Nevertheless, an essential building block of an ultra-reliable wireless system is a model of the wireless channel at the physical layer that captures the statistics of rare events, and extensive fading dips [4], [18].
Current studies on the statistical modeling of the wireless channel operating in the ultra-reliable communication region are categorized into four folds: ) Extrapolation of a wide range of commonly used average statistics-based wireless channel models such as Rayleigh and Rician towards the ultra-reliability region [4], [5]; ) Recommendation of the usage of new channel parameters by incorporating extreme reliability requirements into the communication [22], [18]; ) Usage of a non-parametric statistical learning algorithm to estimate the probability density function (PDF) of the channel distribution [23]; and ) Application of the EVT to characterize the statistics of the channel tail distribution by deriving the statistics of extreme events happening rarely [11], [12].In [4], a simple power-law expression is proposed to extrapolate the cumulative distribution function (CDF) of the commonly used average statistics channels to the ultra-reliable regime of operation.In addition, the authors in [4] apply the power-law results to analyze the performance of receiver diversity schemes and obtain a new simplified expression for MRC applicable for an ultra-reliable communication.In [5], two performance measures, average reliability (AR) for dynamic environments and probably correct reliability (PCR) for the static environments, have been proposed based on the extrapolation-based channels in [4] with the goal of selecting a transmission rate guaranteeing ultra-reliability.However, these extrapolated distributions cannot accurately estimate the lower tail distribution, and therefore, result in several orders of magnitude differences in the estimated PER [11], [12].Also, the channels in [4] are assumed independent and non-identically distributed (i.n.i.d.), which is not a realistic assumption.In [22], [18], the authors have modified the definitions of coherence time and distance to the time and distance over which a channel is predictable with a given reliability, respectively.Accordingly, the authors in [18] use a standard spatially independent quasi-static Rayleigh fading model of wireless channels to examine channel dynamics in the URC context.Nevertheless, none of these frameworks derives ultra-reliability statistics by proposing a channel modeling methodology.In [23], the authors propose a non-parametric statistical learning algorithm that utilizes kernel density estimation (KDE) to estimate the PDF of channel distribution and accordingly, selects a proper transmission channel rate that meets the requirements of ultrareliability.However, the proposed data-driven non-parametric approach requires a massive number of training samples, about 10× −1 , to meet the targeted error probability .Only recently in [11], [12], we have proposed a novel channel modeling methodology based on EVT to derive the statistics of the lower tail distribution while efficiently dealing with a massive amount of data.EVT is a unique and robust framework to develop techniques for modeling the statistics of rare events based on the implementation of mathematical limits as finitelevel approximation [4], [14], [24].
EVT has been used at the data link and network layers to model the tail statistics of queue length and delay [25]- [27] and derive closed-form asymptotic expressions for the throughput, bit error rate (BER), and PER over different fading channels [28]- [31].Moreover, block maxima EVT models have been applied to model the maximum end-to-end (E2E) latency by using the generalized extreme value (GEV) distribution for virtual reality (VR) applications in Terahertz [32].Furthermore, EVT has been used to derive the distribution of the maximum end-to-end SNR in an opportunistic relay selection-based cooperative relaying network consisting of a large number of independent and non-identical relay links [33].Additionally, EVT and federated learning (FL) approaches have been combined in [34] to propose a Lyapunov-based distributed transmit power and resource allocation procedure for vehicular users (VUEs), while the statistics of the queue lengths exceeding a high threshold are characterized by using the GPD.However, these upper layer studies assume the average statistic-based channel models such as Rayleigh, Rician, or Nakagami fading; therefore, their usage may not be suitable in a system operating at URC [4], [11], [12].To fill this gap, in [11], we have used EVT techniques at the physical layer to determine the optimum threshold below which the received power samples are considered extreme events and included in the tail distribution, then model the channel tail distribution by using the generalized Pareto distribution (GPD), and finally, assess the validity of the proposed EVT-based model by using the probability plots.Additionally, in [12], we have extended [11] to model the tail of the non-stationary channel distribution by first determining an external factor causing non-stationarity and, accordingly, dividing the channel data into multiple stationary sequences to apply the EVT techniques initially presented in [11] and model the parameter of the fitting distribution as a change-point function of time.Additionally, recently, in [35], [36], we have proposed a novel EVT-based framework dealing with a relatively low number of data samples to estimate the optimal transmission rate and validate it by assessing the outage probability so that reliability constraints are met with given confidence for ultra-reliable communications.Nevertheless, we have not considered spatial diversity to address the reliability constraints in more than one dimension and assess the dependency between extreme events.Therefore, there is a lack of studies focusing on channel modeling strategies in the ultra-reliable regime of operation by discovering the inter-relation of extreme events.Multivariate extreme value theory (MEVT) is a robust statistical discipline that develops techniques to model the relation of rare events based on the multidimensional limiting relations [4], [14], [24].
The goal of this paper is to propose a novel channel modeling methodology based on MEVT for a system using spatial diversity in MIMO-URC to derive the lower tail statistics of the received signal power in multiple dimensions while efficiently dealing with a massive amount of corresponding data.We focus on the two-dimensional or bi-variate case to highlight the main concepts and issues of MEVT without increasing the complexity of the notation required for an entire multivariate perspective.The modeling approach adapts MEVT to (i) fit GPD to the tail distribution of the received power samples exceeding a given threshold in each sequence and derive the scale and shape parameters of GPD while the thresholds are determined optimally; (ii) apply Fréchet transformation to each data sequence so that the marginal distributions of the resulting variables are Fréchet; (iii) determine the dependency factor between the Fréchet transformed sequences; (iv) apply the logistic distribution approach and Poisson point process approach to fit bi-variate GPD (BGPD) to the tail of the joint distribution of the Fréchet transformed sequences, considering the dependency factor; and (v) assess the validity of the fitted BGPD model by checking the mean constraint on the probability measure function.MEVT extends the ideas of univariate EVT to the analysis of dependent extreme events for multiple random variables by introducing several features of MEVT that are not available in EVT analysis, including modeling dependence structures and modeling the joint distribution of multiple extreme events.The original contributions of the paper are listed as follows: • We propose a comprehensive channel modeling methodology for a system operating at URC based on MEVT for the first time in the literature.The methodology consists of techniques for deriving the lower tail statistics of each channel data sequence by using Uni-variate GPD (UGPD), fitting BGPD to the tail of the joint probability distribution by using the logistic distribution-based and Poisson point process-based approaches, and assessing the validity of these two proposed models by incorporating the probability measure function of the Pichanks coordinates.
• We propose novel techniques based on the logistic distribution to fit BGPD to the tail of a joint probability distribution of channel data, for the first time in the literature.We first determine the dependency factor among the Fréchet transformed sequences, and then derive a closedform expression for the BGPD model.These techniques are original and contribute to the advancement of statistical models for analyzing joint probability distributions of extreme channel samples.• We propose novel techniques based on Poisson point process approach to represent the Fréchet transformed channel tail data using their pseudo-polar Pickands pairs, radial and angular components, for the first time in the literature.We use these representations to determine the probability measure function of the Pickands angular component and the corresponding BGPD model.• We introduce a novel approach for assessing the validity of the derived BGPD models for the joint distribution of the channel tail based on the verification of the mean constraint on the corresponding probability measure function of the Pickands coordinates obtained from BGPD approaches, for the first time in the literature.• Based on the data collected within the engine compartment of Fiat Linea using one transmitter and two receivers under various engine vibrations and driving scenarios, we demonstrate the superiority of the proposed methodology for deriving the tail statistics of multivariate extremes compared to the conventional extrapolation-based models of the average statistics channels to the ultra-reliable region, in terms of the modeling accuracy.
The rest of the paper is organized as follows.Section II describes the system model and assumptions considered through-out the paper.Section III describes the basics of uni-variate and multivariate EVT together with the theorems used in the development of the proposed multivariate channel modeling approach.Section IV presents the proposed channel modeling framework based on MEVT for characterizing the joint multichannel tail distribution in the ultra-reliable region.Section V provides the channel measurement setup and the performance evaluation of the proposed algorithm in determining the optimum threshold, fitting BGPD to the joint probability distribution of the extremes, and comparing the proposed methodology to the conventional methods in terms of the estimation accuracy.Finally, concluding remarks and future works are given in Section VI.

II. SYSTEM MODEL
We consider MIMO for a single transmitter (Tx)-receiver (Rx) pair communicating over a stationary channel, i.e., the parameters of the channel distribution are fixed over a vast period.If the channel is non-stationary according to the Augmented Dickey-Fuller (ADF) test results [12], the external factors causing time variation are determined such that the sequence is divided into groups, each of which can be considered stationary, as explained in detail in [12].Then, the tail distribution of each stationary sequence is modeled by using the GPD.The GPD is used in EVT to estimate the tail of distribution by modeling the probabilistic distribution of the values exceeding a given threshold [14], [24], [37]- [42].Please note that the same procedure can be applied to multiple transmitter and receiver pairs.The transmit power is assumed to be fixed and known in advance.Therefore, estimating the received signal power is equivalent to estimating the squared amplitude of the channel state information [5], [11].
Ultra-reliability is defined as communication with target packet error probability in the range of 10 −9 -10 −5 .We assume that the outage is the only source of packet error and is defined as the received power being less than a predefined threshold [4], [43], [44].Since the main focus is on modeling the channel behavior in the ultra-reliable region by employing the distributions, modeling and validation techniques adopted from extreme value theory in an offline manner, the delay required for collecting a large number of empirical samples is not considered in the first step of the study.The design of a real-time URLLC system based on these statistics requires the consideration of a limited amount of data requiring the inclusion of confidence intervals in the parameter estimation or usage of transfer learning techniques, but it is out of the scope of this paper and subject to future work.
In order to model the multivariate extremes of received powers, the sequences of measured received power samples at Tx-Rx pairs are converted into a sequence of independent and identical distributed (i.i.d.) samples by removing their dependency via declustering approach [11].Upon applying EVT to the resulting sequence of i.i.d.samples in each pair, the optimum threshold is determined for each sample sequence.Determining an optimum threshold is of paramount importance as it determines the number of samples that are included in the channel tail and considered as an extreme value.
Then, the parameters of the Pareto distribution associated with the optimum threshold are estimated by using the maximum likelihood estimator (MLE).A multivariate version of GPD, i.e., a family with which to approximate a joint distribution on regions where we observe joint extreme values, is obtained by using the MEVT.

III. BACKGROUND
EVT is a robust framework that models the probabilistic distribution of extreme events occurring rarely.EVT has been applied in different fields, including hydrology to quantify risks of extreme floods, rainfalls, and waves [24], and finance to estimate losses due to extreme events [37].However, it has been recently used in the telecommunication field to analyze the behavior of extreme values either in network traffic, worst-case delay, queue lengths, throughput, and BER/PER of URLLC [35], [25], [45], or in channel modeling to statistically derive the lower tail of the received power for URC systems [11]- [12].In the following, Section III-A presents the uni-variate technique for modeling the extremes of a single process.Section III-B gives the general modeling technique of multivariate extremes to describe the extreme value interrelationships of two or more processes.Sections III-B1 and III-B2 present bi-variate extreme value techniques based on the logistic family distribution and Poisson point process approach, respectively.Finally, Section III-B3 provides techniques to assess the validity of the proposed bi-variate extreme value modeling techniques.

A. Uni-variate Extreme Value Theory
Uni-variate EVT (UEVT), in general, focuses on the representation and modeling techniques for the extremes of a single process.UEVT is used for modeling extreme events in two main ways.The first concerns models for block maxima by using the (GEV) distribution given by where , , and are the location, shape, and scale parameters of the GEV distribution, respectively [14,Theorem 3.3].The second uses EVT to characterize the tail of a distribution, i.e., the extremes of a sequence, by modeling the probabilistic distribution of values exceeding a given threshold through the generalized Pareto distribution (GPD).Assume that { 1 , ..., } is an i.i.d.stationary sequence, where denotes the ℎ received power for ∈ {1, ..., }, and is the total number of received power samples.Then, according to the EVT, the tail distribution of the power sequence, i.e., the probabilistic distribution of the power values exceeding a given threshold , can be expressed as where is a non-negative value denoting the exceedance below threshold , i.e., ( = − ), denotes any below threshold ; ( ) expresses the GPD; and and ˜ = + ( − ) are shape and scale parameters of the GPD, respectively.It should be noted that and are the location and scale parameters of the GEV distribution fitted to the CDF of = { 1 , ..., }, respectively [11, Theorem 1], [46].

B. Bi-variate Extreme Value Theory
When studying the extremes of two or more processes, the individual process can be modeled by using uni-variate techniques.However, the possible dependency between the extreme events requires the investigation of their joint behavior.Bi-variate EVT (BEVT) allows us to estimate the probability of exceeding thresholds simultaneously and analyze the interdependence of two variables in the extreme value region [47].
In the following, we provide the Theorems and Corollaries required to develop our channel modeling methodology.Theorem 1 incorporates BEVT applications to investigate the general form of BGPD models that are valid to estimate the bi-variate tail distribution.Then, based on Theorem 1, we obtain the logistic family distribution-based BGPD and Poisson point process-based BGPD models in Theorems 2 and 3, respectively.Also, we define the Fréchet transformation in Definition 1 as the input of BEVT in Theorem 1 is required to be Fréchet distributed.Additionally, the definition of Pseudopolar Pickands coordinate transformation is provided in Definition 2 to be used in the Poisson point process approach in Theorem 3.Moreover, the constraints on the Pickands coordinates are defined in Definition 3, which will be used to estimate the probability density function and the probability measure function of the BGPD model in Theorem 2. Finally, Theorems 4 and 5 express the requirements of a valid BGPD model based on the mean constraint on the probability measure functions of the Pickands coordinates corresponding to the BGPD models obtained in Theorems 2 and 3, respectively.Definition 1 (Fréchet transformation).Suppose ( 1 , 1 ), ..., ( , ) are independent realizations of a random variable ( , ) with joint distribution ( , ).For optimum thresholds and , each marginal distribution of has an approximation of the form (2), with the parameters ( ˜ , ) and ( ˜ , ) for and , respectively.Let's apply Fréchet transformation to variables and to induce new variables ˜ and ˜ , given by and where = ( < ) and = ( < ).Then, the joint distribution function of ˜ and ˜ , ˜ , has margins that are approximately standard Fréchet distribution, i.e., where ˜ and ˜ are any realizations from ˜ and ˜ , respectively [14].
, where * ˜ , = max =1,..., { ˜ }/ , and * ˜ , = max =1,..., { ˜ }/ , and ( ˜ , ˜ ) are independent vectors with standard Fréchet marginal distributions as expressed in (3) and (4), respectively, and is the number of realizations in the tail of or .Then, where − → denotes limit of distribution, and is a bi-variate non-degenerate distribution function with the form and is a distribution function on [0, 1] satisfying the mean constraint Any family of distributions that arise as limits in ( 5) can be considered the class of bi-variate extreme value distributions [14], [48].
Proof.Please refer to [49] for the proof.
Theorem 1 implies that the class of bi-variate extreme value distributions is in one-to-one correspondence with the set of distribution functions on [0, 1] satisfying (8).If the probability measure is differentiable with probability density ℎ, integral ( 7) is simplified to This assumes that the thresholds and are small enough to justify the limit (5) as an approximation [14].
Theorem 1 also implies that for generating a class of BGPD models expressed in (6), it is required to obtain a parametric family for over [0, 1] whose mean is 0.5 for every value of the parameter.However, in practice, it is hard to find parametric families whose mean is parameter-free and for which the integral (9) is tractable.Two approaches that can be utilized to model the bi-variate extreme value distribution are the logistic family and the Poisson point process, which will be discussed in the following.
The variables ˜ and ˜ of the logistic distribution family in (10) are exchangeable and have correlation = 1− 2 [24].The value very close to 1 denotes strong dependence between the variables and , even at moderately extreme levels.
2) BGPD Based on Poisson Point Process: Definition 2 (Pickands coordinates).Let ( ˜ , ˜ ) be independent vectors with standard Fréchet marginal distributions as expressed in (3) and ( 4).Let's define the transformation ( ˜ , ˜ ) by where ∈ , is the length of vector ˜ or ˜ .( ˜ , ˜ ) is called transformation with respect to the Pickands coordinates ( , ), where is pseudo-polar angular component measuring angle in [0, 1] scale, and is pseudo-polar radial component measuring the distance from the origin [14].This transformation is one-to-one with the inverse The pseudo-polar Pickands mapping has the same geometrical interpretation as standard polar coordinates with the difference that polar coordinates use the Euclidian norm for the angular and radial component, while the Pickands coordinates use the sum norm [49]- [50].
Proof.Let and denote the Pickands coordinates of ˜ and ˜ given by where ˜ and ˜ have standard marginal Fréchet distributions.Then, the intensity function of on space in Theorem 3 is [14], [47], [50] ( , Additionally, let be the point process defined as (12) with intensity function ( , ) in (16).Then, Λ( ) in Theorem 3 is given by where is the maximum of − ˜ / and − ˜ / 1− , depending on how is defined based on Remark 1.The Pickands probability measure ( ) has the following property [50]- [53]: The fact expressed in Remark 1 will be used in Theorem 5 to assess the validity of the Poisson point process-based BGPD.
3) BGPD Model Assessment: Any distribution function defined on the space [0, 1] in (7) that satisfies the mean constraint in (8), gives rise to a valid limit in (5) [14].Therefore, if ( ˜ , ˜ ) is a valid model to estimate the tail of bivariate extremes, its corresponding probability measure function ( ) should satisfy the 0.5 mean constraint according to Theorem 1.

2) Conditional on
> 0 , the radial component of the Pickands coordinates is uniformly distributed.Therefore, Theorem 4. If the BGPD model based on the logistic distribution family ( ˜ , ˜ ), given by (10), is a valid model to estimate the tail of the bi-variate extremes, the mean constraint defined in (8) should be satisfied on ( ) = ∫ ℎ ( ) , where = ∫ ( )d > 0, and ( ) is the Pickands density function given by where is the BGPD expressed as (10) and 0 is the optimum cut-off point obtained based on the constraints on the Pickands coordinates in Definition 3.
Proof.Please refer to [49] and [50] for the proof.

Theorem 5. If the BGPD model
( ˜ , ˜ ) based on the Poisson point process approach given by ( 13) is a valid model to estimate the tail of the bi-variate extremes, the mean constraint defined in ( 8) is always satisfied for the ( ) function in (14).
Proof.Let and denote the Pickands coordinates of ˜ and ˜ based on Definition 2.Then, referring to (18), we have where ∫ [0,1] ( ) = 1, and therefore, the 0.5 mean constraint defined in ( 8) is satisfied for the probability measure function ( ).

IV. PROPOSED CHANNEL MODELING METHODOLOGY
We propose a novel BEVT-based channel modeling methodology with the goal of estimating the joint lower tail statistics of multiple channels in the ultra-reliable regime.The methodology consists of the following steps: The sequences of measured received power samples are converted into sequences of i.i.d.samples by removing their dependency via declustering, where the samples are divided into multiple clusters, each of which includes consecutive dependent observations, and clusters are separated by a specific sample gap to ensure the independency between clusters [11].Upon applying EVT to the resulting sequence of i.i.d.samples, the optimum thresholds are determined.Since the bi-variate analysis is restricted to the time intervals in which both thresholds are exceeded, we keep only the exceedances that occur in the same time interval.Otherwise, we ignore the exceedance in any sequence.Thereafter, the parameters of the UGPD associated with the optimum thresholds are estimated by using the MLE.Then, we assess the validity of the fitted UGPD model to the tail distribution by using the probability plots, including the probability/probability (PP) plot and the quantile/quantile (QQ) plot.Afterward, we model the inter-relationship of bi-variate extremes based on Theorem 1 by applying two methods, logistic distribution and Poisson point process.In the logistic distribution approach, upon determining the dependency factor between the Fréchet transformed variables in Definition 1, we model the tail of the joint PDF based on Theorem 2 and verify the model according to Theorem 4 by using the doubled-transformed data to the Pickands coordinate based on Definition 2. In the Poisson point process approach, based on Theorem 3, we transform the data sequence in two steps: First, based on the Fréchet transformation in Definition 1, and second, based on Pickands coordinates in Definition 2.Then, we determine the probability measure of the Pickands angular coordinate ( ) and the intensity measure of Poisson point process Λ( ).Finally, we assess the validity of the Poisson point process-based bi-variate model based on Theorem 5.The proposed algorithm for the bi-variate extremes is depicted in Fig. 1 and explained in detail next.A. Modeling the Uni-variate Extremes 1) Declustering Approach: In the declustering approach, first, we assume a low enough threshold and start the first cluster by the first sample below this threshold.Then, all the successive samples below threshold followed by the first sample of the initial cluster are assigned to the first cluster.Right after observing a sample over , we let the cluster continue for more consecutive values and then, if no value below is detected, close the cluster.The next cluster starts with the following value below the threshold .Upon determination of the clusters of samples, we extract the minimum value of each cluster, apply EVT to the i.i.d.cluster minima, and model their tail distribution by using GPD.The declustering approach is based on the fact that the minimum values of the clusters are far enough to be considered independent and identically distributed [11], [14].
2) Optimum Threshold Determination: Determination of the optimum threshold is of paramount importance as it specifies where the tail of the distribution starts and distinguishes between extreme events happening rarely and nonextreme values.The optimum threshold of each uni-variate case is determined by using the mean residual life (MRL) and parameter stability methods.Based on the MRL method, is the optimum threshold if is the highest threshold below which the expected value of samples exceeding is a linear function of , i.e., ( − | < ) is linear against .The remarkable advantage of the MRL method is its simplicity, as it can be applied to the data sequence prior to the estimation of the UGPD parameters.However, due to the lower precision of the MRL method compared to the parameter stability method, it is sometimes difficult to obtain the optimum threshold explicitly.Therefore, the MRL method is usually utilized as a complementary method or in the case that the optimum threshold determination does not require high precision.The parameter stability method states that the optimum threshold is the highest threshold below which the estimated shape and modified scale parameters of the UGPD fitted to the tail distribution associated with a variety of thresholds is a linear function of the threshold.It is worth noting that the modified scale parameter is defined as * = ˜ − .Additionally, the linearity relation is assessed by using the R-squared statistical measure, denoted by 2 [11]- [12].
3) Model Validity Assessment: The validity of the GPD model is assessed by using probability plots, i.e., probability/probability (PP) and quantile/quantile (QQ) plots.These plots are graphical techniques used to assess the validity of the models fitted to the empirical values.In the PP plot, we plot the empirical CDF of the occurrence of an extreme value versus the corresponding CDF obtained by the GPD, while in the QQ plot, we plot the empirical extreme quantile versus the corresponding value obtained by the inverse of GPD [11].If the GPD appropriately models the extreme values exceeding threshold , then both PP and QQ plots should fit the unit diagonal line, i.e., the 45 • line [14], [37].

B. Evaluation of the Correlation Coefficient
The first step in the bi-variate extreme analysis is to assess the amount of dependence in the tails [47].However, since the bi-variate analysis is restricted to those time intervals in which both thresholds and are exceeded, before determining the correlation among the tail samples, we are required to revise the tail samples by removing those exceedances happening in one sequence but not in the other one, at the same time interval [14].Accordingly, we consider a time interval consisting of samples and check if and are exceeded simultaneously in sequences and , respectively.If so, we capture the minima of clusters; otherwise, we ignore the minima of each cluster [14], [47].Then, EVT is applied to the obtained exceedances to determine the UGPD parameters.
Upon determining the extremes of each sample sequence, if the extremes of the two sequences are independent, no bivariate modeling like ( 5) is required.Otherwise, the existence of correlation among the tail samples suggests the requirement of investigating the inter-relationship of the extreme values of receivers.On the other hand, correlation among the received powers in the total samples is used to assess the feasibility of the spatial diversity.If the correlation coefficient among the total samples is between 0.1 and 0.5, the spatial diversity is suggested [54].Otherwise, the employment of spatial diversity is pointless due to the concurrent fading at different links.If the spatial diversity is feasible and the correlation among the tail samples is high enough, we are required to analyze the bi-variate tail characteristics for a URC system.

C. Modeling the Multivariate Extremes Based on the Logistic Distribution Approach 1) Data Conversion Based on the Fréchet Transformation:
Applying the Fréchet transformation by using the Definition 1 on the obtained tail sequences and from Section IV-A, we induce variables ˜ and ˜ whose marginal distribution functions have Fréchet distribution for < and < , approximately.This transformation is required as the input of the bi-variate extreme modeling approach is expected to have the Fréchet distribution [14], [24].
2) Determining the Dependency Factor : To determine the dependency parameter of the logistic model in ( 11), we either estimate the correlation coefficient between the variables ˜ and ˜ and then, calculate as 1 − , or directly estimate it by MLE, where the log-likelihood function is defined as where, and ( ˜ , ˜ ) is the logistic distribution function defined in (11).
3) BGPD Model Based on the Logistic Distribution: The bi-logistic family of distributions formulated in ( 11) is used to model the bi-variate extremes based on Theorem 1. Upon obtaining the dependency parameter , as well as the Fréchet transformed variables ˜ and ˜ , we build function (11) and then, by taking the exponential of its negative function, we determine the BGPD as expressed in (10).If this BGPD model is reliable to characterize the tail of the bi-variate distribution, the mean constraint on the probability measure function defined in Theorem 4 should be satisfied.

D. Modeling the Multivariate Extremes Based on the Poisson
Point Process Approach 1) Pseudo-polar Pickands Coordinates Transformation: In Poisson point process-based approach for modeling the bivariate extremes, it is required first to transform the data from Cartesian to the polar coordinate: ( ˜ , ˜ ) → ( , ).To this end, we apply pseudo-polar Pickands transformation based on Definition 2 that induces two new variables, radial component and angular component .
2) Determining the Probability Measure Function of the Radial and Angular Components: The CDF of the probability measure function of the angular component ( ) is determined based on Theorem 3.This probability measure will be used in the next step to determine the intensity of Poisson point process Λ( ) and, later, to assess the validity of the proposed Poisson point process-based BGPD, according to Theorem 5.
3) BGPD Model Based on the Poisson Point Process: Upon determining the angular component probability measure function ( ), we compute the density function of the Poisson point process Λ( ) for the defined space based on Theorem 3. Λ( ) is actually equivalent to ( ˜ , ˜ ) in the logistic distribution based-BGPD model.Afterward, the BGPD model based on the Poisson point process approach is determined as exp(−Λ( )).

E. BGPD Model Validation
1) Model Validation for Logistic-based BGPD: According to Theorem 4, if the logistic-based BGPD is a valid model to estimate the tail of the bi-variate extremes, the probability measure function ( ), conditional on { : > 0 }, needs to satisfy the 0.5 mean constraint, i.e., ∫ ( ) = 0.5, where is the Pickands angular component of ˜ and ˜ .
To determine 0 based on Definition 3, a plot of versus for < 0 , referred to as the -plot, can be utilized to address the first constraint in Definition 3 by checking the dependency between and .The dependency of and is assessed based on the correlation results in which and are independent if their correlation is less than the critical value 0.05.Additionally, to address the second constraint in Definition 3, the distribution of the radial components 0 < < 0 is expected to be fitted to the uniform distribution, where ), where and are the number of transmitters and receivers, respectively, and is the number of the training samples for individual channel sequences.

V. NUMERICAL RESULTS
The goal of this section is to evaluate the performance of the proposed channel modeling algorithm in estimating the BGPD model fitted to the joint distribution of the channel data sequences based on two approaches: the logistic distribution approach and the Poisson point process approach, and also compare their performances with the traditional extrapolationbased approaches to estimate the statistics of the channel tail for a system operating in the spatial diversity in MIMO-URC.In the traditional extrapolation-based approach, upon estimating the distribution of the existing channel data for the reliability order of 10 −3 -10 0 PER [55]- [56] for individual channel data sequences, we determine their joint probability distribution, and then, extrapolate it towards the ultra-reliable region 10 −9 -10 −5 PER [4].
We have collected channel measurement data within the engine compartment of Fiat Linea under various engine and driving scenarios at 60 Gigahertz (GHz) by using a Vector Network Analyzer (VNA) (R & S ® ZVA67).The VNA is connected to the transmitter and receivers through the R & S ® ZV-Z196 and PE361 port cables, respectively, with 610 millimeter (mm) length as shown in Fig 2 .The output power range at the port can be decreased to −100 dBm to measure deep fading.The transmitter is an omnidirectional antenna operating from 58 GHz to 63 GHz with 0 decibel isotropic (dBi) nominal gain.The receivers are horn antennas operating between 50-75 GHz with a nominal 24 dBi gain and 11 • and 9.5 • horizontal and vertical half power beamwidth, respectively.The antennas are connected to the coax cables through the waveguide to the coax adaptor operating at the frequency span of 50 −65 GHz, with insertion loss 0.5 decibel (dB) and impedance 50 Ohm (Ω).
The data were collected while driving the car at the Koc University campus and emulating different scenarios, including starting/stopping the car, moving up/down on a ramp, and driving on a flat road.More than 10 6 successive samples are captured from each receiver antenna for about 1.5 hours with a time resolution of 3 ms.We use MATLAB for analyzing the data and the implementation of the proposed algorithm as well as the traditional extrapolation-based approach [4].Meanwhile, the estimation error of the proposed methodology and the extrapolated-based approach have been reported by means of the Root Mean Square Error (RMSE) metric.
In the following, first, in Section V-A, we provide the numerical results in the determination of the optimum threshold over which the tail statistics are derived for two univariate channel data sequences obtained from receivers Rx1 and Rx2, and then, validate the tail model corresponding to the optimum threshold by means of the probability plots.Upon determining the correlation among the tail samples in Section V-B, we apply the logistic distribution approach to model the tail distribution of the bi-variate extremes in Section V-C.Next, in Section V-D, we estimate the tail distribution of the bi-variate extremes by applying the Poisson point process approach.Afterward, in Section V-E, we assess the validity of the determined BGPD models based on logistic distribution and Poisson point process approaches.Finally, we compare the performance of both proposed methodologies, logistic distribution-based and Poisson point process-based approaches, to that of the traditional extrapolation-based technique in the estimation of the bi-variate channel tail statistics in Section V-F.

A. UGPD Model
Towards determining the optimum threshold of the GPD fitted to the tail distribution of samples in each group, the declustering approach is used to remove the dependency among the samples and obtain i.i.d.observation for the EVT input.In this regard, the MRL plot and parameter stability method are utilized to determine the thresholds and , and the minimum gaps(mgs) between the samples mg 1 and mg 2 , for receiver Rx1 and Rx2, respectively.We skip the illustration of these results and only present the critical information as the subject of the uni-variate channel tail modeling has been discussed in our previous paper [11] in detail.According to the results, for the uni-variate channel data sequence obtained from receiver Rx1, = −15 dBm is the optimum threshold below which the 2 > 0.95 for all mg 1 > 1.Additionally, = −30 dBm is the optimum threshold below which the linearity condition is observed for all mg 2 > 1 for the univariate channel data of receiver Rx2.Upon determining the optimum thresholds ( and ) and minimum gaps mg 1 and mg 2 between the clusters, the probability plots, including the PP plot and QQ plot, are used to validate the accuracy of the UGPD models fitted to the channel tail distribution of the i.i.d.samples obtained from receivers Rx1 and Rx2, corresponding to the optimum thresholds and , respectively.

B. Determination of Correlation
To determine the correlation among the tail samples, first, we have considered the time intervals consisting of 1000 samples, and then, at each interval, obtained those minima of Rx1 and Rx2 sequences that simultaneously exceed = −15 dBm and = −30 dBm, respectively.The parameters of UGPD fitted to these new exceedances are ( , ˜ ) = (−0.1469,4.0367), and ( , ˜ ) = (−0.4245,8.5886) for receivers Rx1 and Rx2, respectively.The correlation among the total samples obtained from Rx1 and Rx2 is about 0.2766, which confirms the applicability of spatial diversity as it is less than 0.5.Additionally, the 0.3957 correlation coefficient among the tail samples indicates the existence of correlation in the tail samples of Rx1 and Rx2 and, therefore, confirms the necessity of studying the inter-relationships of the extreme values.

C. BGPD Based on the Logistic Distribution Approach
Fig. 3 illustrates the CDF of the transformed data by using the Fréchet transformation for the received power tail samples obtained from receivers Rx1 and Rx2.For small values of ˜ and ˜ , the Fréchet transformation is not able to fit appropriately to the empirical data.However, as the transformed variables ˜ and ˜ increase and correspond to the extreme values with extreme fading, the Frećhet distribution fits well to the empirical values.This confirms that ˜ and ˜ are the reliable transformed variables for building the BGPD model.Fig. 4 shows the proposed BGPD model determined based on the logistic distribution approach, ( ˜ , ˜ ).Empirical joint CDF of the tail samples has also been depicted for better compatibility assessment of the logistic-based BGPD.The BGPD model obtained by the logistic distribution approach can estimate the empirical joint CDF with RMSE 0.8655.Additionally, the estimation accuracy increases as we approach more extreme values, which confirms the ability of the proposed methodology to model the worse extreme values.However, the exact pattern of the empirical CDF has not been preserved by the logistic-based BGPD model.

D. BGPD Based on the Poisson Point Process Approach
Fig. 4 also depicts the proposed BGPD model determined based on the Poisson point process approach.It can be seen that the estimated joint CDF based on the Poisson point process approach is in good agreement with the corresponding empirical CDF with an RMSE of 0.8686.Although the Poisson point process-based approach for modeling the bi-variate tail distribution seems to have the same performance as the logistic distribution-based approach according to the RMSE results, it preserves the shape of the empirical joint CDF very well.This is mainly due to the higher complexity of the

E. Model Validity Assessment
Figs. 5 and 6 are utilized to determine the optimum cut-off point 0 for which the Pickands constraints are satisfied and ( ) is defined on > 0 .Fig. 5 illustrates the -plot of the Poisson point process approach for determining the optimum cut-off point 0 .In this figure, the radial component of the Pickands coordinate is plotted versus the corresponding angular component, where and are the pseudo-polar Pickands coordinates of ˜ and ˜ .Please note that {( , )| < 0 } correspond to independent pairs and {( , )| > 0 } are matched to dependent pairs of and .According to this plot, 0 ≈ −0.47 is the minimum radial component above which the correlation coefficient of Pickands radial and angular components is less than the critical value 0.05, resulting in independent and pairs.Therefore, 0 ≈ −0.47 is the cutoff radial component that satisfies the first Pickands constraint.Fig. 6 shows the distribution of the radial component of the Pickands coordinate for > 0 where is the pseudopolar Pickands radial coordinate of ˜ and ˜ , and 0 ≈ −0.47 is the optimum radial cut-off point captured from Fig. 5.It is observed that the distribution of the radial components > 0 , i.e., those radial components independent from their angular component pairs, fits the uniform distribution, which means that ( > ) = 0 , and results in satisfying the second constraint on Pickands coordinates.
We assess the validity of the proposed BGPD models shown in Fig. 4 based on the CDF results of Pickands probability density functions ℎ ( ) and ℎ ( ) obtained from the logistic-based for > 0 , and Poisson point process-based BGPD approaches, respectively.Accordingly, it is observed that the CDF of the estimated density function of the angular Pickands coordinate , for ℎ ( ), the 0.5 mean constraint in (8)  based BGPD approach, the 0.5 mean constraint is satisfied as ∫ ( ) = 0.5064 for the derived ( ).

F. Comparison with Conventional Tail Estimation Method
Fig. 4 additionally illustrates the CDF of the normalized power values for the traditional extrapolation-based approach compared to the empirical and the proposed BGPD models based on the logistic distribution and Poisson point process approaches.To derive the traditional extrapolation-based model, we first extract the observations from Rx1 and Rx2 corresponding to the reliability order of 10 −3 -10 0 PER.Then, upon fitting a variety of distributions to the samples in each observation sequence, the best-fitting distribution is determined as the Gaussian distribution according to the Akaike information criterion/Bayesian information criterion (AIC/BIC) metric.It is worth noting that Gaussian and Rician distribution had almost the same AIC/BIC for both sequences.However, for the sake of simplicity and without loss of generality, we opt to continue with the Gaussian distribution in our analysis as its joint PDF calculation is more straightforward.Then, we determine the joint bi-variate CDF of the Gaussian distribution for variables and , corresponding to the received power values obtained from Rx1 and Rx2, respectively.Upon determining the joint CDF of and , we extrapolate it towards the ultrareliable region based on the tail approximation approach in [4].According to the results illustrated in Fig. 4, the proposed BGPD models outperform the extrapolated Gaussian model remarkably for a URC system, where the reliability order is in the range of 10 −9 -10 −5 .The RMSE of the extrapolated bivariate Gaussian model is 7.8177 × 10 09 , indicating the superiority of the proposed model with a significant improvement in the estimation of the tail distribution for the ultra-reliability region, where the reliability orders are in the range of 10 −9 -10 −5 .

VI. CONCLUSIONS
In this paper, we introduce a novel framework based on the extreme value theory with the goal of estimating the multivariate channel tail statistics for a MIMO-URC.The proposed methodology utilizes MEVT to model the tail of the joint probability distribution by using two methods for modeling bi-variate extreme values: the logistic distribution approach, which uses the logistic distribution to determine dependency factors and obtain a model, and the Poisson point process approach, which estimates the probability measure function of the Pickands angular component to model the bi-variate extreme values.The proposed multi-dimensional channel modeling methodology achieves a significantly better fit to the empirical data in the lower tail than the conventional extrapolation-based approach.In the future, we are planning to extend this work by determining the estimated BGPD parameters and optimum transmission rate for a delayconstrained real-time communication system leading to the design of a MIMO-URLLC system.Since the derivation of the statistics requires a large amount of data, such real-time system design requires either the inclusion of confidence intervals in the parameter estimation or the adoption of transfer learning techniques, such as knowledge-assisted training, where the digital twin of a real network is built with network topology, channel and queueing models for offline training, which is then fine-tuned in the real environment with a smaller amount of data.
Model Validation for Poisson Point Process-based BGPD: According to Theorem 5, the Poisson point processbased BGPD is, by default, a valid model to estimate the tail of the bi-variate extremes as its intensity function ( ) satisfies the mean constraint.Therefore, the constraint ∫ ( ) = 0.5 is insured in Poisson point process-based BGPD while determining ( ).The complexity of the proposed MEVT-based channel modeling framework is (

Fig. 2 .
Fig. 2. Transmitter and receiver antennas in the engine compartment.

Fig. 3 .
Fig. 3. CDF of the received power tail samples transformed based on the Fréchet transformation: (a) for the transformed samples of receiver Rx1, → ˜ , and (b) for the transformed samples of receiver Rx2, → ˜ .

Fig. 5 .Fig. 6 .
Figs. 5 and 6 are utilized to determine the optimum cut-off point 0 for which the Pickands constraints are satisfied and ( ) is defined on > 0 .Fig.5illustrates the -plot of the Poisson point process approach for determining the optimum cut-off point 0 .In this figure, the radial component of the Pickands coordinate is plotted versus the corresponding angular component, where and are the pseudo-polar Pickands coordinates of ˜ and ˜ .Please note that {( , )| < 0 } correspond to independent pairs and {( , )| > 0 } are matched to dependent pairs of and .According to this plot, 0 ≈ −0.47 is the minimum radial component above which the correlation coefficient of Pickands radial and angular components is less than the critical value 0.05, resulting in independent and pairs.Therefore, 0 ≈ −0.47 is the cutoff radial component that satisfies the first Pickands constraint.Fig.6shows the distribution of the radial component of the Pickands coordinate for > 0 where is the pseudopolar Pickands radial coordinate of ˜ and ˜ , and 0 ≈ −0.47 is the optimum radial cut-off point captured from Fig.5.It is observed that the distribution of the radial components > 0 , i.e., those radial components independent from their angular component pairs, fits the uniform distribution, which means that ( > ) = 0 , and results in satisfying the second constraint on Pickands coordinates.We assess the validity of the proposed BGPD models shown in Fig.4based on the CDF results of Pickands probability density functions ℎ ( ) and ℎ ( ) obtained from the logistic-based for > 0 , and Poisson point process-based BGPD approaches, respectively.Accordingly, it is observed that the CDF of the estimated density function of the angular Pickands coordinate , for ℎ ( ), the 0.5 mean constraint in (8) is satisfied as ∫ ( ) = 0.5001, indicating the validity of the BGPD model based on the logistic distribution approach.On the other hand, referring to the Pickands density function ℎ ( ) obtained from the Poisson point process-

Niloofar
Mehrnia received a B.S. degree in biomedical engineering and M.S. and Ph.D. degree in electrical and electronics engineering from Isfahan University, Istanbul Sehir University, and Koc University in 2014, 2017, and 2022 respectively.From 2019, she also worked as a researcher at Koc University Ford Otosan Automotive Technologies Laboratory (KUFOTAL).Her research interests include wireless channel modeling, ultra-reliable and lowlatency communications (URLLC), extreme value theory (EVT), and vehicular communications.Sinem Coleri is a Professor and the Chair of the Department of Electrical and Electronics Engineering at Koc University.She is also the founding director of Wireless Networks Laboratory (WNL) and director of Ford Otosan Automotive Technologies Laboratory.Sinem Coleri received the BS degree in electrical and electronics engineering from Bilkent University in 2000, the M.S. and Ph.D. degrees in electrical engineering and computer sciences from University of California Berkeley in 2002 and 2005.She worked as a research scientist in Wireless Sensor Networks Berkeley Lab under sponsorship of Pirelli and Telecom Italia from 2006 to 2009.Since September 2009, she has been a faculty member in the department of Electrical and Electronics Engineering at Koc University.Her research interests are in 6G wireless communications and networking, machine learning for wireless networks, machine-to-machine communications, wireless networked control systems and vehicular networks.She has received numerous awards and recognitions, including N2Women: Stars in Computer Networking and Communications, TUBITAK (The Scientific and Technological Research Council of Turkey) Incentive Award and IEEE Vehicular Technology Society Neal Shepherd Memorial Best Propagation Paper Award.Dr. Coleri has been Interim Editor-in-Chief of IEEE Open Journal of the Communications Society since 2023, Executive Editor of IEEE Communications Letters since 2023, Editor-at-Large of IEEE Transactions on Communications since 2023, Senior Editor of IEEE Access since 2022 and Editor of IEEE Transactions on Machine Learning in Communications and Networking since 2022.Dr. Coleri is an IEEE Fellow and IEEE ComSoc Distinguished Lecturer.