{'fraction_non_alphanumeric': 0.0358870026348264, 'fraction_numerical': 0.015075209285261649, 'mean_word_length': 4.619190404797601, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Integration of max-stable processes and Bayesian model averaging to predict extreme climatic events in multi-model ensembles

Projections of changes in extreme climate are sometimes predicted by using multi-model ensemble methods such as Bayesian model averaging (BMA) embedded with the generalized extreme value (GEV) distribution. BMA is a popular method for combining the forecasts of individual simulation models by weighted averaging and characterizing the uncertainty induced by simulating the model structure. This method is referred to as the GEV-embedded BMA. It is, however, based on a point-wise analysis of extreme events, which means it overlooks the spatial dependency between nearby grid cells. Instead of a point-wise model, a spatial extreme model such as the max-stable process (MSP) is often employed to improve precision by considering spatial dependency. We propose an approach that integrates the MSP into BMA, which is referred to as the MSP-BMA herein. The superiority of the proposed method over the GEV-embedded BMA is demonstrated by using extreme rainfall intensity data on the Korean peninsula from Coupled Model Intercomparison Project Phase 5 (CMIP5) multi-models. The reanalysis data called APHRODITE (Asian Precipitation Highly-Resolved Observational Data Integration Towards Evaluation, v1101) and 17 CMIP5 models are examined for 10 grid boxes in Korea. In this example, the MSP-BMA achieves a variance reduction over the GEV-embedded BMA. The bias inflation by MSP-BMA over the GEV-embedded BMA is also discussed. A by-product technical advantage of the MSP-BMA is that tedious `regridding' is not required before and after the analysis while it should be done for the GEV-embedded BMA.


Introduction
The multi-model ensemble has the potential to combine the strengths of individual models and to better represent forecast uncertainty than the use of a single model. Various studies have shown the advantage of merging forecasts from multiple climate models for climate projection (e.g., Casey 1995; Coelho et al. 2004;Stephenson et al. 2005; Kharin et al. 2013; Arisido et al. 2017). A commonly used multi-model ensemble method is Bayesian model averaging (BMA), which assumes a probability distribution around the forecasts of individual models to establish the prior and uses a calibration period to determine the static weights for each individual model (Parrish et al., 2012;Wang et al., 2017).
When our interest is in forecasting extreme climatic events, the generalized extreme value distribution (GEVD) is typically used as an assumed probability distribution in BMA. The GEVD encompasses all three possible asymptotic extreme value distributions predicted by large sample theory (see, e.g., Leadbetter et al. 1983). It is widely used to analyze univariate extreme values (Coles, 2001). Zhu et al. (2013) performed a GEV-embedded BMA of extreme rainfall intensities by using seven regional climate models (RCMs). Hereafter, this method is referred to as the GEV-BMA. In that study, the method was applied to each grid cell separately; hence, it ignored the spatial dependency between nearby grid cells. Consequently, the GEV-BMA may be less reliable for small samples and could suffer from local fluctuations, which may be heavily influenced by outlier observations. Thus, it may be better to generalize this separate GEVD approach to a model taking into account the spatial dependency in the BMA framework. Therefore we consider a spatial extreme model, namely a max-stable process (see, e.g., de Haan 1984;Smith 1990; ). Max-stable process (MSP) is an infinite dimensional extension of the GEVD. The approach considered here, an MSP-embedded BMA, is a spatial extension of the GEV-BMA proposed by Zhu et al.(2013). Hereafter, the proposed method is referred to as the MSP-BMA.
MSPs are widely used for modeling spatial extremes, which deal with the spatial dependency of extreme events arising between nearby locations. The MSP estimates not only the GEVD parameters as continuous functions of space but also spatial dependent structure. The advantage of the MSP is that it can be used to address problems concerning the aggregation of global processes over many locations. Because it builds a model for the whole set of grid cells, it avoids local fluctuations that may be heavily influenced by outlier observations. Thus, more precise and robust return levels of extreme events can be obtained.
The present study aims to determine whether the integration of the MSP and BMA can add value when predicting future extreme values using multiple simulation models. It was motivated by Parrish et al. (2012) and Madadgar and Moradkhani (2014), who integrated data assimilation methods and copulas into BMA, respectively. Their results demonstrated that the predictive distributions after these integrations are more accurate and reliable, less biased, and more confident (i.e., having lower uncertainty). Sabourin et al. (2013) considered using BMA for multivariate extreme events. More recently, Najafi and Moradkhani (2015) and Yan and Moradkhani (2016) considered spatial hierarchical Bayesian modeling of extremes to obtain spatially distributed BMA weights corresponding to each simulation model. Our approach is in the same line of integrating spatial models into BMA framework. We are unaware of any work in which the MSP has been integrated into BMA. This is a reason why we considered MSP-embedded-BMA approach in this study. Based on the success of their research, we also expect a similarly good result by integrating the MSP and BMA.
To demonstrate the prediction ability of our approach, we use the annual maximum daily precipitation as simulated by Coupled Model Intercomparison Project Phase 5 (CMIP5) models in historical and future experiments over the Korean peninsula. Based on the variance estimates of return levels and other criteria, our approach seems outperforming the GEV-BMA.
The remainder of this paper is organized as follows. Section 2 describes models of extreme events including the GEVD and the MSP. Section 3 explains the GEV-embedded BMA and introduces the integrated MSP-BMA. Section 4 provides a numerical study of extreme rainfall as simulated by CMIP5 models over the Korean peninsula. A computationally faster version of the MSP-BMA is considered in Section 5. Section 6 concludes with discussion.

Extreme value distribution
The GEVD is widely used to analyse univariate extreme values. The three types of extreme value distributions are sub-classes of GEVD. The cumulative distribution function of the GEVD is as follows: where µ, σ, and ξ are the location, scale, and shape parameters, respectively. The particular case for ξ = 0 in Eq. (1) is the Gumbel distribution, whereas the cases for ξ > 0 and ξ < 0 are known as the Fréchet and the negative Weibull distributions, respectively ( [6]). Assuming the data, the annual maxima of daily precipitation in this study, follow (approximately) a GEV distribution, the parameters can be estimated by using the maximum likelihood method (see, e.g., Wilks 1995 or Coles 2001) or the method of L-moments estimation (Hosking, 1990). The maximum likelihood estimator is less efficient than the L-moments estimator in small samples for typical shape parameter values (Hosking et al. 1985). The L-moments method is used in this study because relatively short 30-year samples are analyzed for each comparison period. Moreover, the formulae used to obtain the L-moments estimator are simple compared with obtaining the maximum likelihood estimator, which needs an iterative computation until convergence. Those areξ = 7.8590c + 2.9554c 2 where and l 2 and l 3 are the sample second and third L-moments (Hosking, 1990). The other parameters are then given by where l 1 is the first sample L-moment. It can be helpful to describe changes in extremes in terms of changes in extreme quantiles. These are obtained by inverting (1): where G(z p ) = 1 − p. Here, z p is known as the return level associated with the return period 1/p, since the level z p is expected to be exceeded on average once every 1/p years (Coles, 2001). For example, a 20-year (50-year) return level is computed as the 95th (98th) quantile of the fitted GEVD. This return level is used as the main quantity in the numerical study to compare the GEV-BMA with the MSP-BMA in the subsequent sections.

Models for spatial extremes
Extreme value distributions are typically applied independently for the datasets of each weather station (grid cell, in this study). However, a better approach may be to build a model for all the weather stations together to avoid local fluctuations and reflect the spatial dependency between observations from nearby stations. A variety of statistical tools have been used for the spatial modeling of extremes including Bayesian hierarchical models, copulas and max-stable processes. The first model is often referred to as "latent variable" approach (Davision et al., 2012). Spatial hierarchical models introduce spatial variation in the parameters of the extremal response distributions. By letting the parameters of the response distributions depend on the unobserved latent processes, the spatial hierarchical model induces dependence in responses by integration over the latent process. A Gaussian process is frequently used for the latent process.
This approach is common in geostatistics with nonnormal response variables (Diggle and Ribeiro, 2007), and because of the complexity of the integration involved is most naturally performed in a Bayesian setting, using Markov chain Monte Carlo algorithms to perform inferences (see, e.g., Sang and Gelfand, 2010;Cooley et al., 2012;Fawcett and Walshaw, 2016). Compared to the models based on copulas and MSP, it is reported that the spatial hierarchical modeling makes quantifying uncertainty more straightforward, and it allows a better fit to marginal distribution, but it fits the joint distributions of extremes poorly (Davision et al., 2012).
Copula provides a unifying framework to modeling multivariate data, and so may be useful to model spatial extremes. Given marginal GEVDs and some correlation functions, various copulas including Gaussian and Student t are available. Some copulas satisfying the maxstability property, called as extremal copulas, may be useful to model spatial extremes (see, e.g., The max-stable processes (MSP) is an useful natural extension of the extremal (univariate and multivariate GEV) models. It is an infinite-dimensional generalization of the univariate GEV model, and is particularly applicable in the context of a time series or a spatial process. This theory was introduced by De Haan (1984) and has been extended by a number of researchers such as Smith (1990) . It, however, also bring challenges. For example, the likelihood is hard to tractable for even moderate sample size problems, and consequently composite (or pairwise) likelihood methods are utilized.
For future projections and uncertainty assessment of extreme events from an ensemble of climate models, many proposed the integration of spatial extreme models and BMA framework. We are, however, unaware of any work in which the MSP has been integrated into BMA. Thus we propose an integration of MSP and BMA in multi-model ensemble, in this study.

Max-stable processes
The max-stable process (MSP) Z(·) is defined as the limit process of the maxima of the independent and identically distributed random fields Y i (x), x ∈ R d (De Haan, 1984) and can be expressed as if a limit exists for all x in R d with normalizing constants a n (x) > 0 and b n (x) ∈ R. For a fixed point (x 1 , · · · , x d ) in space, any d-dimensional marginal distribution of the MSP belongs to a class of multivariate extreme value distributions. In practice, this means that the resulting parameters µ(x), σ(x) > 0, and ξ(x) are the continuous functions to be estimated. Thus, the MSP naturally permits modeling and predictions with data-level spatial dependence. We then build the regression-based forms for the parameters µ(x), σ(x), ξ(x) and estimate the spatial dependence parameters ( [48]). These regression-based forms, however, may introduce some estimation bias for MSP while it are convenient to predict something for non-observed sites. The precision of the estimated parameters in the MSP is expected to be improved over the point-wise (or at-site) GEV model because the MSP allows more data to be used in the model fitting by combining data from multiple sites across a spatial domain (Ribatet, 2013;Gaume et al., 2013;and Oesting and Stein, 2017).
Given a series of n observations at K spatial locations, the aim of a statistical analysis would be to fit an MSP by using assumed forms for the parameters µ(x), σ(x), and ξ(x), while also estimating spatial dependence. However, for more than K = 2 spatial locations, the distribution function of the general MSP has no analytically tractable form, which thereby presents a problem for practical statistical model fitting. Therefore, ad hoc characterization methods have been proposed. Smith (1990) introduced the first, and Schlather (2002) and Kabluchko et al. (2009) later introduced other methods. Schlather's characterization involves a correlation function ρ(x), deriving a bivariate cumulative distribution function as where h is the distance between x 1 and x 2 . Further, ρ ∈ [0, 1], where the lower bound corresponds to the independence of extremes. To take into account the spatial dependency, we used a powered exponential correlation function ρ(h) of distance h, in this study [8]: where τ and η are the scale and shape parameters that control the range and smoothness of the process, respectively. The Brown-Resnick model [20] uses the following variogram [8]: The full likelihood of the MSP is often analytically unknown or computationally prohibitive to calculate. The following pair-wise log likelihood function is used because only the bivariate density function is specified; where x i(k) is the k-th observation at the i-th site for i ∈ {1, · · · , d}. The use of the pair-wise likelihood method leads to an underspecified model. Thus, instead of the Akaike information criterion, the Takeuchi information criterion (TIC) is used to select a better model for the estimated parametersψ = (μ(x),σ(x),ξ(x)) [44]: where tr(A) stands for the trace of matrix A, and where ∇ means the derivative with respect to ψ. A model with a low TIC is a good one. In standard BMA, the conditional distribution of each individual model is assumed to follow a normal distribution, which is valid for temperature and sea-level pressure, for example. For other variables such as maximum precipitation or extreme wind speed, the GEVD might be a better alternative for representing the distribution of the model's output. In their climate change study of extreme rainfall, Zhu et al. (2013) used the GEVD in a BMA framework, where the weight of each forecast model was calculated by comparing the reanalysis data with the historical data from a simulation model. Moreover they applied the bootstrap technique instead of employing the expectation maximization algorithm to compute the weights of the models and variance of the BMA predictions of rainfall intensity. Now, we describe the details of their method, which is referred to as the GEV-BMA in this paper.
Given the relatively short record of 30 years from the climate model and reanalysis data in this study, the parameter estimation of the GEVD may be imprecise and problematic for predicting events beyond the range of the observed data (e.g., return level for 100 years). Thus, Zhu et al. (2013) proposed a bootstrap-based approach to establish the uncertainty of the prediction. Based on the ability of each climate model's historical runs in matching the reanalysis data for each location, the weights for all RCMs were calculated. For each RCM and each bootstrap realization, the GEV parameters were generated and rainfall intensity (i.e., return level) in return period T estimated. The variance in rainfall intensity was calculated for all bootstrap realizations (from i=1 to B). They then used the Gaussian likelihood function as follows: where I i (T ) is the intensity of the i-th bootstrap from the reanalysis data, I k i (T ) is the intensity of the i-th bootstrap from the historical data of model M k , and is the variance in intensity based on the reanalysis data, whereĪ(T ) = 1 . Note that the data from all models were downscaled to have the same grids as the reanalysis data in advance of fitting the GEVD at each site.
If I = I(T ) is the rainfall intensity of return period T predicted by a set of K alternative models, then its distribution conditioned on the reanalysis data R is (Hoeting et al., 1999) where p(M k |R) is the predictive probability of I for model M k and p(M k |R) is the posterior probability of M k . One way of estimating the posterior model probability is to use Bayes' theorem: where p(R|M k ) is the likelihood of model M k and p(M k ) is the prior probability of M k . The above quantity (17) is used as the weights of the BMA as in the below equation (19). Here, p(M k ) may be assigned based on prior information from expert elicitation, professional assessments of the incorporated models, and so on. If such information is unavailable, the prior probabilities could be set to be equal for all the models. In the numerical study of this article, equal priors are used for all 17 models. The likelihood of model M k is approximated by using the generalized likelihood uncertainty estimation method (Beven, 2006;Rojas et al., 2008): The mean and variance of the BMA prediction of rainfall intensity I over K models are then given by where w k = p(M k |R) is the weight for each model k. For the historic data, these quantities are estimated by the calculation from the B bootstrap samples used in (14): and For the future simulation data, the above equations (21) and (22) are not used, but the intensity estimates obtained by fitting GEVDs to each model's future data are used. That is, where the subscript f stands for the future data. For estimation of V ar(I|R, M k ), we may need to consider other quantities such as asymptotic variance.

Bayesian model averaging using MSP
In Zhu et al. (2013), the dataset from a simulation model in each grid point was treated independently even though nearby observations were aggregated. That is, the GEVD was fitted at-site and the return levels in period T were estimated for each grid. The weighted average of these return levels over a suite of simulation models is the GEV-BMA prediction. As mentioned in subsection 2.3, the precision of these estimated return levels can be improved by applying the MSP instead of using point-wise GEV models. Now, we propose an integration of BMA and the MSP with the expectation of an improved ensemble prediction compared with the GEV-BMA. We refer to this new integrated approach as the MSP-BMA.
In the MSP-BMA, the MSPs are fitted to the reanalysis data and the historical data of model M k , respectively, over the spatial domain. Then, the return level, I(T ), at each site x is obtained by using the following equation: where p = 1/T . This is a spatial extension of equation (6). Here, the GEV parameters µ(x), σ(x), ξ(x) obtained from the MSP fitting are functions of the spatial variables given in (26) in the next section. To compute (14), bootstrap samples are used as in the GEV-BMA. These are obtained by random sampling with replacement over the spatial domain, from the reanalysis data and the historical data of model M k . It was sampled 'year-wise' over spatial stations. If the year 1980 is selected, for example, then the observation vector of all stations on the year is used as a bootstrap realization. Then, the BMA weights are calculated by using the same formula as in (17). Moreover, the BMA prediction and uncertainty estimates over multi-models for the future data are calculated by the same equations as in (19) and (20), respectively. The main advantage of the MSP-BMA is that the precision of the BMA prediction is improved, as shown by the illustrative numerical study in the next section. As an advantage of MSP-BMA as a by-product, we have experienced a technical convenience of handling data in applying MSP-BMA to multiple climate models compared to applying GEV-BMA. It is because that each climate model has different resolution and so has different locations and sizes of grid cells. Moreover, those are also different from those of reanalysis data. Thus, to compare each model to reanalysis data and to combine multi-models, we needed a pre-processing to make the center points of grid cells to be same (i.e., regridding) across multi-models and reanalysis data in applying GEV-BMA, while we do not need it in applying MSP-BMA.

Comparison methods
To compare the performances of the GEV-BMA and the MSP-BMA, we considered the following two criteria. The first one is the quantile per quantile (qq) plot for the group-wise maxima to check the goodness-of-fit (Davison and Gholamrezaee, 2012). Here, the empirical groupwise maxima for each year are selected from the given observations across the sites, while the theoretical group-wise maxima are obtained from the average of the B randomly generated MSPs for each year. This is done for all the years (30 years in this study) to construct a maxima series. Then, the quantiles from each maxima series are plotted. This plot with a confidence band is drawn for each method and compared.
The second criterion is the variance in return levels of each method for each site. This is calculated by using equation (22) for both methods. Lower variance values represent less uncertainty in the estimation. This comparison is conducted for each site (10 locations in this study). The relative improvement is calculated by

Numerical study 4.1 Datasets
We analyze annual maximum daily precipitation (AMP1) as simulated by the CMIP5 models in the historical experiments  and the experiments for the 21st century (2021-2050). We considered the models only employing one radiative forcing scenario, called Representative Concentration Pathways (RCP) 4.5. Table 1 lists the model names, institutes, and resolution of the 17 models considered in this study. For the reanalysis data, we use APHRODITE (Asian Precipitation Highly-Resolved Observational Data Integration Towards Evaluation, v1101) daily gridded precipitation data. This is the long-term (1951 onward) continental-scale daily product that contains a dense network of daily rain-gauge data for Asia including the Himalayas, South and Southeast Asia, and mountainous areas in the Middle East (Yatagai et al., 2009). The grid resolution of APHRODITE is 0.5 o × 0.5 o for East Asia including Korea. We treat this reanalysis data as the "ground truth." In a statistical analysis such as the GEV-BMA, we need to match the grid point of the reanalysis data to that of the CMIP5 models. However, the resolution of the CMIP5 models varies, as seen in Table 1, and thus we need a standard grid to efficiently use these datasets. For such a standard grid, we set the grids as 2 o × 2 o , as seen in Figure 1. For each box in the figure, there are 16 observations of AMP1 from APHRODITE. We treat the average of these 16 observations as the representative AMP1 of the grid box. For our study, 10 boxes are chosen over the Korean peninsula. The red triangle in Figure 1 stands for the center site of the box. For each site, a time series of 30 years (from 1971 to 2000) AMP1 is used. For the historical and future data from each simulation model, the Kriging method is applied to obtain the observation (AMP1) at the center site. The "DiceKriging" package in R software (Roustant et al., 2012) was used with Matern 5/2 covariance structure. Now, we have three kinds of datasets for each site: reanalysis and historical data from 1971 to 2000 (period 0) and future simulation data from 2021 to 2050 (period 1).

MSP modeling
For the data analysis, we employ a linear transformation to make the station's geographical information be in [0, 1]:  To compute the MSP, we use the R package "SpatialExtremes" developed by Ribatet et al. (2013). From the reanalysis data, we build regression-based forms to the location and scale parameters by including latitude and longitude, where x denotes the station. We set ξ(x) as a constant ξ, which is typical in MSP modeling (Buishand et al., 2008;Westra and Sisson, 2011;Lee et al., 2013), because we are studying with a small number of grid cells. The selected model at the minimum of the TIC is as follows: where lat(x) and long(x) are the latitude and longitude of the site x and µ 0 , µ t , µ g , µ t2 , µ tg , σ 0 , σ g , σ tg are the regression coefficients corresponding to the variables. This selected regression form is fixed for the MSP for the historical and future data in this study because the model selection based on the TIC for the 17 CMIP5 models took huge computing time. However, the regression coefficients are estimated for every dataset by maximizing the pair-wise log likelihood function given in (11). Here it is notable, as commented by a reviewer, that the polynomial form for GEV parameters is just a convenience assumption but may not be very flexible. Treating GEV parameters as Gaussian processes in spatial Bayesian hierarchical model may provide a more flexible alternative. To select the characterization among the three approaches (Smith, Schlather, and Brown-Resnick), we check the extremal coefficient plots for each simulation model. Figure 2 shows typical plots drawn for some of the simulation models: CNRM-CM5, NorESM1-M, IPSL-CM5A-LR, and MRI-CGCM3. From these plots, we choose the Brown-Resnick characterization in this study. Figure 3 shows the group-wise maxima qq plots. The left panel is from the GEV fitting and the right one is from the MSP fitting to the reanalysis data. The right plot provides a better qq than the left one in the sense that the points in the right plot are closer to the straight line and more located inside the confidence band. This means that the MSP is more acceptable than the GEVD for this dataset.

BMA result for comparison
For each site, the variance estimates of the 20-year return levels are calculated from the GEV-BMA and the MSP-BMA by using the reanalysis and historical data. Figure 4 compares these two kinds of variances (unit: mm 2 ). The Y-axes and X-axes stand for the variances estimated by using the GEV-BMA and the MSP-BMA, respectively. The GEV-BMA has greater variance than the MSP-BMA for eight of the 10 sites. This may mean that the MSP-BMA provides more precise predictions with less uncertainty than the GEV-BMA. Table 2 shows the estimated variances of the 20-year return levels obtained by using both BMA methods. The relative improvement is calculated by the formula (24) in which the MSP-BMA achieves a 42% relative variance reduction on average over the GEV-BMA. It may not be surprise that the variances obtained from MSP is smaller than point-wise GEV, but it could also due to the spatial modeling bias and the variance estimation may be too small to be true. Thus we used a bootstrap method to estimate the bias more reasonably. Bootstrap samples are obtained by random sampling with replacement over the spatial domain, from the reanalysis data and the historical data of model M k . It was sampled 'year-wise' over spatial stations, as we did in the previous section 3.2. The bias is calculated by the difference between the true 20-year return level and estimated levels obtained by a bootstrap sample. This calculation is repeated B times, and the average of B differences is considered to be the bias of the model. Here, for the true 20-year return level, we used the estimate of 20-year return level obtained by applying GEV-BMA model for the original data. If we set the true 20-year return level differently from this study, the bias estimates for both models may be changed . The  table 3 shows the bias estimates calculated by using GEV-BMA and MSP-BMA methods. The table shows that the MSP-BMA has bias inflation over GEV-BMA at eight sites. The bias of MSP-BMA may be introduced by assuming a regression model for each GEV parameter. By considering both results from Tables 2 and 3, we know that it is an example of bias-variance trade-off problem, as commented by a reviewer. We leave more study on this matter to future work.

Summary and discussion
In this study, a MSP-embedded BMA is proposed to project extreme climates based on multimodel ensemble simulations. The proposed method, referred to as the MSP-BMA, is a spatial extension of BMA embedded with a point-wise GEVD. The superiority of the proposed method over the GEV-embedded BMA is demonstrated by using extreme rainfall intensity data on the Korean peninsula from CMIP5 multi-models. The APHRODITE reanalysis data and 17 CMIP5 models are examined for 10 grid boxes over the Korean peninsula. A conclusion of this study is that the MSP-BMA outperforms the GEV-embedded BMA by generating lower prediction variance estimates for extreme rainfall intensity data on the Korean peninsula. The MSP-BMA achieves a 42% relative variance reduction over the GEV-embedded BMA on average. However, the regression-based forms in MSP may introduce some bias, which needs more study on this bias-variance trade-off problem. A by-product technical advantage of the MSP-BMA is that tedious 'regridding' is not required before and after the analysis while it should be done for the GEV-embedded BMA.
As commented by a reviewer, the constant ξ parameter assumption in MSP effectively reduces the estimation uncertainty while there is no such a constraint from point-wise GEV method. For a fair comparison, one can also modify the point-wise GEV method to a constant shape parameter GEV model (Wang et al. 2016). We leave this matter to future work.