On the Smoothing of the Generalized Extreme Value Distribution Parameters Using Penalized Maximum Likelihood: A Case Study on UVB Radiation Maxima in the Mexico City Metropolitan Area

: This paper concerns the use and implementation of penalized maximum likelihood procedures to ﬁtting smoothing functions of the generalized extreme value distribution parameters to analyze spatial extreme values of ultraviolet B (UVB) radiation across the Mexico City metropolitan area in the period 2000–2018. The model was ﬁtted using a ﬂexible semi-parametric approach and the parameters were estimated by the penalized maximum likelihood (PML) method. In order to investigate the performance of the model as well as the estimation method in the analysis of complex nonlinear trends for UVB radiation maxima, a simulation study was conducted. The results of the simulation study showed that penalized maximum likelihood yields better regularization to the model than the maximum likelihood estimates. We estimated return levels of extreme UVB radiation events through a nonstationary extreme value model using measurements of ozone (O 3 ), nitrogen oxides (NO x ), particles of 10 µ m or less in diameter (PM 10 ), carbon monoxide (CO), relative humidity (RH) and sulfur dioxide (SO 2 ). The deviance statistics indicated that the nonstationary generalized extreme value (GEV) model adjusted was statistically better compared to the stationary model. The estimated smoothing functions of the location parameter of the GEV distribution on the spatial plane for different periods of time reveal the existence of well-deﬁned trends in the maxima. In the temporal plane, a presence of temporal cyclic components oscillating over a weak linear component with a negative slope is noticed, while in the spatial plane, a weak nonlinear local trend is present on a plane with a positive slope towards the west, covering the entire study area. An explicit spatial estimate of the 25-year return period revealed that the more extreme risk levels are located in the western region of the study area. decrease in the intensity of UVB radiation in the central region, which return to the initial levels in 2018 and 2019. These results show that our model allowed us to identify the complete temporal dynamics of the trend throughout the study period. This is one of the strengths of the proposed model, which


Introduction
Ultraviolet radiation can cause different effects on Earth's life. In living organisms, UVB radiation destroys DNA, produces protein denaturation, triggers coagulation of albumin, as well as erythema and skin problems. In humans, UVB radiation causes the weakening of immune system, creates as well as in environmental, anthropogenic and geophysical processes. Particularly, it has been used to study the long-term risks in rainfall [21], winds [22], heat waves [16] and earthquakes [23,24]. In these studies, it can be seen that the trend of extreme values has been adjusted using several approaches. In fact, one of these approaches is that if the patterns follow the law determined by a model, then the GEV parameters of the corresponding model are estimated [25,26]. In contrast, for some others, it is more appropriate to adjust the trend with smoothing functions [21,27]. In both cases, the trend is adjusted by estimating the parameters on predictors of the location parameter of the GEV distribution. Analogously, an adjustment similar to the logarithm of the scale parameter is made [28][29][30]. The contrast occurs for the shape parameter, which is assumed to be constant because the estimation is numerically fraught when this parameter is allowed to be too flexible [28]. Extreme nonstationary values have also been studied extensively using the Bayesian approach. For instance, Gaetan and Grigoletto [15] studied rainfall maxima with Markov random fields approximated based on smoothing kernel, Reich et al. [16] studied heat waves using a Bayesian hierarchical model with the generalized Pareto distribution (GPD) and Sang and Gelfand [31] studied the extreme values of spatial stochastic processes and modeled the observed trend as a function of spatial covariates.

Study Area
The

A Nonstationary GEV Model
Inferences about the parameters of the extreme values can be made with the exact distribution of the maximum of a random sample when the cumulative distribution function of the target population is known. However, in large samples the exact distribution function tends to concentrate the mass of the probabilities in a single point, known as degenerate distribution. This kind of distribution is not useful for further analysis. In other cases neither the population distribution is known, nor the sample is small. One solution to these limitations is to approximate the asymptotic distribution of the extremes through the limit distribution of a properly rescaled sequence. Consider Y 1 , ..., Y n a sample of independent and identically distributed random variables with cumulative distribution function F Y (y) and let M n = max (Y 1 , ..., Y n ), then the only limiting nondegenerate distribution G n = (M n − a n )/b n as n → ∞ (if such a sequences of constants {b n } and {a n } such that for each n ∈ N, b n > 0 exist) is the generalized extreme value (GEV) distribution [11]: Here, µ ∈ R, σ > 0 and κ ∈ R are the location, scale and shape parameters, respectively. The quantile function of the GEV distribution, obtained with Q (p) = G −1 (p), is given by The GEV distribution was derived using the stationarity assumption inherited from a random sample. In a real scenario, the maxima are usually not identically distributed, i.e., the mean of the distribution varies as a function of covariates. In such cases, we establish a predictor to the parameters of location, scale and shape of the form µ t = µ (X t1 , ..., X tk ), σ t = σ (X t1 , ..., X tk ) and κ t = κ (X t1 , ..., X tk ) [20].

Proposed Approach
Similar to the approach of generalized linear vector models proposed by Yee and Stephenson [28] and the analysis of nonstationary extreme values proposed by Coles [20], we associated a linear predictor to the parameters of the GEV distribution. The linear predictor expresses the relationship of a set of covariates with the maxima through the parameters, which usually consists of linear functions. The structure of the linear predictor is analogous to linear regression models based on spline-based functions. However, the regression is not directly done on the response variable, but is assigned a linear function with a radial basis kernel to approximate the trend of UVB radiation maxima through the location parameter. We chose a flat function for the scale and shape parameter respectively, because the estimation is numerically fraught when this parameters are allowed to be too flexible. Therefore, the proposed model is as follows: where φ (·) is a real value function, σ t , κ and µ t are scale, shape and location parameters, respectively, x t = x t1 , . . . , x t(p 1 −1) is the vector of covariates for the t-th observation, scaled and centered, and k j corresponds to the j-th centroid obtained by the method of hierarchical clustering among x t , t = 1, 2, . . . , n [32]. Note that the set of location parameter can be expressed in matrix notation as µ = Xβ (1) + Zu (1) where β (1) = c 0 , c 1 , . . . , c p 1 −1 is a vector of size p 1 , u (1) = d 1 , . . . , d p 2 is a vector of size p 2 , X is a n × p 1 matrix which has the additional column of one for the coefficient of c 0 and Z ij = exp φ x i − k j 2 is a p 2 × n matrix of kernel basis which consists of a set of columns obtained by using radial Gaussian function used to approximate the trend and capture the interactions between covariates.

Penalized Maximum Likelihood
Let y = (y 1 , ..., y n ) be a sample of n extremes. The maximum likelihood estimator for the nonstationary GEV is defined as the estimator that maximizes the joint density of the n random variables: It can be seen that maximizing the likelihood function with respect of the GEV parameters is equivalent to maximizing the log-likelihood function: is a vector of p × 1 parameters, the linear predictor of the location parameter can be written as µ = Cb (1) . The penalty of the parameters is introduced through the following matrix: where I p 1 and I p 2 are identity matrices of order p 1 and p 2 respectively, σ 2 β , σ 2 u , σ 2 v and σ 2 κ are values that control the degree of shrinkage on coefficient estimates and D d u 1 = ∆ d u 1 consists of the vector of dth differences of u 1 .
Therefore the penalized log-likelihood of the model is: The gradient of the likelihood is given by: where c t = x t , z t . The Hessian matrix is given by H = n ∑ t=1 H (t) as follow: Considering ⊗ as the usual Kronecker product, the elements of the Hessian matrix are: We can observe that the gradient and Hessian equations have been expressed in terms of those derived from the density function of the response variable, in this case the GEV distribution. Thus, we considered this formulation to be advantageous because these expressions could be used analogously in different smoothing models with penalties that have a linear predictor similar to Equation (2).

Simulation Study
In order to examine the performance of the semi-parametric GEV model defined in Equation (2), a simulation study was conducted. The nonlinear system that is addressed by Equation (5) was used to simulate the trend of n = 500 extreme values sampled from a nonstationary GEV distribution. In this function, x 1 and x 2 , correspond to the longitude and latitude, respectively. To perform a more realistic simulation, the intervals for x 1 and x 2 were chosen on the current range of the UVB radiation data. Therefore, the x 1 values were generated with randomly spaced data on the interval [98.88, 99.38], and the values of the covariate x 2 were randomly selected from the interval [19.15, 19.73], with The performance of proposed model was evaluated through the simulation of a simple function which consists of two critical points where the function has a maximum and minimum, respectively. We simulated the data using the inverse transform method proposed by Ross [33] as follows: (1) we generated random values for x 1 and x 2 and obtained µ t , σ and κ, using Equation (5); (2) we generated a random value q with uniform distribution and (3) simulated values were obtained through the inverse of the GEV distribution, given by the Equation (1), using the value of q obtained in step 2 and the values of µ t , σ and κ obtained in step 1. Regarding the number of knots, we consider two settings, the first one using a set of basis functions with p 2 = 20 knots and the second with p 2 = 80 knots, in both cases we choose σ 2 β = 1, σ 2 u = 1, σ 2 v = 1 and σ 2 κ = 1. We can see in Figure 2, that the log-likelihood has stabilized at settings p 2 = 80. Once the model was estimated, the estimates obtained for the shape parameter were −0.4000 and −0.4016, respectively. The estimate of the scale parameter was 0.1000 and 0.1002, respectively. The true functions of µ expressed as a function of the covariates x 1 and x 2 are shown in Figure 3a. The function corresponding to µ involves the sum of the exponentials to the square of the covariates x 1 and x 2 , which represents a nonlinear surface in the spatial plane, similar to the conditions in the extreme values that can occur in real situations of nonstationary extreme values.
As expected, Figure 3c shows that the estimation improved by increasing the number of knots in the spline functions. The function for σ and κ are constant in the covariate space. The goodness of fit for the proposed model using the simulated data was evaluated using the mean square error (RMSE) and Pearson's correlation, given by Equations (6) and (7). We simulated a new data set of 500 observations of model in Equation (2) and then calculated the RMSE with respect to their estimated values. For the location parameter, in the model with k = 80, we obtained a RMSE of 0.0788 and a correlation of 0.9949 between the predicted values and the testing data. Regarding the interpretation of the coefficients, due to the nonlinear form of the estimated function, the estimated coefficients of the spline model cannot be interpreted directly as marginal changes of the covariates, however, we can still analyze the adjusted function in the covariable space and obtain enough information about the behavior of the trend of extreme values for different values of the covariates.
where µ t is the t-th simulated trend,μ t is its corresponding estimate;μ t andμ t are the means of the set of simulated and estimated values, respectively; and σ µ t and σμ t are their corresponding standard deviations. An inspection of the estimated functions presented in Figure 3 show that the proposed model with both p 2 = 20 and p 2 = 80 recovers the original form of the real function used in the simulation. In the case of the estimated smoothing function build with 20 knots, i.e., p 2 = 20, we observed a more pronounced border effect, similar to the case of spline functions when the response variable has a normal distribution. This border effect is visibly reduced as we increased the number of nodes, in the model with k = 80 of Figure 3c, in which we graphically observed that the shape of the estimated function is more similar to the true function Figure 3a. Similarly, a comparison of estimators of σ in both p 2 = 20 and p 2 = 80 reveals that the estimation of the scale parameter of the extreme values based on the model given in Equation (2) improves in the case when the number of nodes is increased when k = 80. In contrast, the estimators of the shape parameter are slightly skewed to the right in both cases.

Data Description
The data corresponded to 397 observations of bi-monthly maxima of UVB radiation and its corresponding atmospheric covariates, obtained between 1 January 2000, and 30 September 2018, at 7 fixed monitoring stations of the Red Automática de Monitoreo Atmosférico, RAMA (http://www.aire. cdmx.gob.mx/default.php). This monitoring network subsystem is one of the three subsystems of the Sistema de Monitoreo Atmosférico (SIMAT) established by the Comision Ambiental Metropolitana of Mexico City to monitor compliance with ambient air quality standards. The information obtained by the measuring instruments of the RAMA network is concentrated in a computer that sends the information continuously through the modem to the Control Center. The gases are measured in real time, by different methods. O 3 is measured by photometry in the ultraviolet range; NO x by chemiluminescence; CO by nondispersive spectroscopy by correlation; relative humidity is measured using a sensor, capacitor-type, of polymer thin film.

Data Analysis
The extreme values were obtained from space-time using the block maxima methodology. The statistical analysis was performed using the R 3.4.2 software [34]. In the spatial plane, we considered each station as a block, a total of 7 stations were used. In the temporal plane, the width of the time interval was two months. Due to problems in the measurement instruments, we identified 195 records which had missing data in any of the covariates associated with a maximum of radiation.
These observations were excluded from the study. We estimated a GEV model for the maxima of UVB radiation on the Mexico City metropolitan area using a multivariate smoothing functions with spatio-temporal and environmental covariates like latitude (s 1 ), longitude (s 2 ), time (t), ozone (O 3 ), nitrogen oxides (NO X ), carbon monoxide (CO), relative humidity (RH), particulate matter of 10 µm or less in diameter are called PM 10 (PM10) and sulfur dioxide(SO 2 ), grouped into the X matrix to fit the trends in the nonstationary GEV model. In order to avoid abrupt changes between the coefficients, we assign a P = 1 σ 2 1 I P 2 + D d D d penalty matrix to the coefficients u 1 s in the model in Equation (4), where D d (with d = 1) is a matrix such that D d u 1 = ∆ d u 1 constructs the vector of dth differences of u 1 .

Results and Discussion
The descriptive statistics of the data by each monitoring station are shown in Table 1. On this table, we can see that there are differences in the distributions of the extreme values in each of the monitoring stations, i.e., the mean of the distributions is not constant. We verified these results in the boxplot presented in Figure 4. By a simple inspection of the descriptive statistics, we observed that the distribution of the extremes is not stationary, therefore the use of a nonstationary model of extreme values is justified for the analysis of trends. The inspection of Table 1 indicates that in station 1, located in the Acatlán area, the maximum intensity recorded was 6.09 W/m 2 , in contrast to the maximum intensity measured in station 5, located in the San Agustín area, located in the municipality of Ecatepec de Morelos, in the State of Mexico, where the maxima was 5.65 W/m 2 . Similarly, the station 3, located at the Merced and the station 7, located at Tlalnepantla, showed the lowest UVB radiation values in comparison with the other stations. Moreover, on these three stations where we observe more frequently intense periods of air pollution within the study region, the level of UVB radiation is lower than other less polluted areas in the MCMA. The results of the comparison between the modeling of maximum likelihood method and the modeling of penalized maximum likelihood is shown in Table 2. In this table, we observed that the estimated parameters have a considerable shrinkage, which is a desirable characteristic that indicates a strong regularization of the model. We validated the results obtained by observing that the value of −12.49 corresponding to the log-likelihood of the adjusted model has been significantly improved in relation to the value of −180.07 corresponding to the log likelihood of the stationary model. We also validated the proposed model using the Deviance statistic [20]. In the contrast of two models, M 1 with θ 1 a parametric vector against another model M 0 with θ 0 a subset vector such as M 0 ⊂ M 1 , the deviance statistic defined by D = 2 (l * n (M 1 ) − l * n (M 0 )), where l * n (M) is the maximized log likelihood function of model M, is used to prove the superiority of the M 1 model. Values of D greater than the quantile 1 − α of the χ 2 distribution with k degrees of freedom, are considered significant, where k is the difference between the dimensions of M 1 and M 0 . In Table 2 shows the results of the validation of the proposed model using the statistical deviance test, which evidences the improvement of the proposed model compared to the stationary model. These results indicate that the deviance statistics is significant with a reliability level of 99%. Therefore, we concluded that our model allows us to explain spatial and temporal trends using the relationship between covariates and the amount of UVB radiation measured at ground surface.  The estimates obtained for the σ and κ parameters were 0.2504 and −0.0356, respectively. The spatial smoothing for the years 2000, 2005, 2010, 2015, 2018 and 2019 for the location functions are shown in Figure 5. The results show well-defined patterns related to the trend in the spatial plane. The magnitude of the UVB radiation maxima decreased as we moved toward the east direction of the study area. This region coincides with the most industrialized areas of the MCMA. Therefore, these results indicate that the air pollution covariates reduce the net amount of UVB radiation that reaches the ground surface. In contrast, we observed that in the less polluted areas of MCMA there is a greater amount of UVB radiation. However, we were able to observe that, although the general trends are maintained throughout the study period, between 2010 and 2015 there was a small decrease in the intensity of UVB radiation in the central region, which return to the initial levels in 2018 and 2019. These results show that our model allowed us to identify the complete temporal dynamics of the trend throughout the study period. This is one of the strengths of the proposed model, which allows us to identify patterns in the distribution of maximum UVB radiation, make inferences and obtain conclusions about extreme values throughout the study region over time. The results also show the advantage of using a spline model with radial-based functions to estimate trends in extreme values. The nonlinear spatial function estimates in each of the different periods through a single model show the existence of the spatial variation of UVB radiation maxima. The proposed model also has the advantage that it includes the effects of covariate interactions in the model through the use of spline functions that depend on the norm of Euclidean distance between covariates and knots. This feature of the model combined with the penalty of the parameters results in a smooth continuous surface as shown in Figure 5. Future research could include the study of other types of distances between observations. The results of the temporal trend of UVB radiation over the years keeping the other covariates constant are shown in Figure 6. In this figure, we can observe the existence of cyclic temporal patterns in the trends of the maxima in the regions located around the monitoring stations. An important finding related to the temporal behavior of the extremes, which can also be seen in Figure 6, is the decrease in the location parameter over time. An explanation for this finding could be the increase in pollution, specifically the amount of ozone (O 3 ), nitrogen oxides (NO x ), particles of 10 µm or less in diameter (PM 10 ) and carbon monoxide (CO), which decreases the amount of UVB radiation as a result of direct chemical reactions or by radiation blockage. The spatial distribution of the extreme values of UVB radiation is influenced by the physicochemical interactions it has with covariates such as ozone, nitrogen oxides, particles of 10 µm or less in diameter (PM 10 ), carbon monoxide (CO), relative humidity (RH) and sulfur dioxide (SO 2 ) [6]. The atmospheric concentrations of some of these covariates also present seasonal behaviors which modify the intensities of UVB radiation over time. These covariates are used in the nonstationary extreme value model to estimate the trend on UVB radiation maxima, through the linear predictors corresponding to the location and scale parameters. In order to increase the likelihood of the model, we built the design matrix by using a nonlinear function of the square of distance of each observation to knots on the vector space of the observations. Each knot represented one of the k centroids resulting from a hierarchical clustering. There are several approaches to obtain a basis in the column space of the covariates, however, considering the sample size and the number of nodes, the radial basis functions are sufficient to obtain a linearly independent set.
One of the most important applications of the models obtained with the analysis of nonstationary extreme values, consists in the elaboration of risk probability maps and the return level maps. The return level Z p is the threshold at which an extreme value is exceeded with probability p, which is expected to occur once every 1/p years (Fawcett and Green [35]). Figure 7 shows the maximum expected UVB radiation for a return period of 25 years. Isolines on the map (Figure 7) were used to visualize the spatial risk of maximum UVB radiation. In fact, the highest UVB radiation values over a 25-year return period can be expected in the west part of the study area in the regions surrounding the SFE and PED monitoring stations (Figure 7). An interesting fact that explains the spatial trend of the maxima is the amount of atmospheric pollution resulted from the emissions of internal combustion vehicles and industrial emissions, among others. Further, in densely populated areas, such as the Merced or Hangares, where a large number of vehicles circulate daily, as well as in industrialized areas such as Naucalpan and Tlalnepantla, we can expect to have the lowest return levels of UVB radiation of the entire study area. An opposite situation occurs in regions farther from urban areas, in which an increase in the estimates of UVB radiation intensity can be observed. Therefore, the map of return levels shown in Figure 7 allows us to confirm the potential risk of dangerous levels of UVB radiation in the study region. These findings should encourage the creation of policies and the revision of standards related to the protection against UVB radiation as well as the delimitation of critical areas of risk. We agree with the results of Ailliot et al. [36], who reported that for stationary case, the imposing constraints improves the performance of the estimation on κ parameter. However, we verified these results in the nonstationary case. Similar to the results of Martins and Stedinger [19], we obtained implausibly large estimates of κ for unconstrained maximum likelihood. We observed that the Newton-Raphson method does not reach the global optimal solution for most of the initial values. This happens because we have enough variables to estimate and the quadratic approximation, which is the basis of the optimization algorithm, is not appropriate to approximate the log likelihood function when the initial values are distant from the optimal value. On this case, the optimization of the penalized log likelihood was carried out in two stages. The first stage consists of finding an initial point or seed, which will be used in the second stage of the algorithm. To achieve this, the optimization was performed in a smaller parametric space. Once the maximum is found in a parametric space with a dimension smaller than the original, we used these values as seeds or initial values to perform the approximation using the Newton-Raphson algorithm in the initial parametric space. We also agree with the results of Coles and Dixon [37], which found that estimators are improved using the maximum penalized likelihood method by restricting the range of κ.
Similar to the work of Bais et al. [38], we conclude that there is a relationship between UVB radiation and ozone, SO 2 and clouds on the spatial UVB radiation distribution across the metropolitan area of Mexico City. Other researchers have found similar relationships through direct chemical studies, concluding that chemical reactions that involve UVB radiation to produce compounds such as ozone, nitric oxide or sulfur dioxide decrease the amount of UVB radiation that reaches the ground [3,39]. There also exists other factors that interact with UVB radiation. In heavily polluted regions, there are several types of particulate matter which block the UVB radiation path [38]. The goal of our study was to analyze the spatio-temporal distribution of the UVB radiation maxima, since extremes have a strong impact on public health, while in most studies only study continuous measurements. We take advantage of the chemical interactions between air pollution and UVB radiation to model the temporal dynamics on the spatial distribution of maximas in the study area.
The monitoring and analysis of UVB radiation levels is a priority concern in terms of public health for all the largest population centers. In previous studies on UVB radiation in the metropolitan area of Mexico City, Acosta and Evans [40] found that UVB radiation levels, measured over the international standard units, reached dangerous levels for humans. We agree with their findings, in which they also detected a strong attenuation of UVB radiation at ground level in the urban troposphere under polluted conditions. However, in contrast to them, we have used the GEV distribution. An alternative to this distribution is the skew generalized extreme value distribution (SGEV) [41], which showed that it improves the return level estimation in the case of a slow convergence or in the heavy-tailed case. Future work should consider the use of the SGEV distribution for analysis of UVB radiation extremes.

Conclusions
In this study, we have developed a nonstationary extreme value model for UVB radiation maxima on the metropolitan area of Mexico City using a semi-parameterized model to obtain a spatio-temporal smoothing of the location parameter of the GEV distribution. We have estimated return levels of extreme events of UVB radiation through a nonstationary extreme value model in which we use both spatial and environmental covariates. UVB maxima were obtained in each of the monitoring stations through the block maxima method. The spatial and temporal trend was approximated by means of the location parameter of the GEV distribution using linear predictors based on Gaussian basis functions of the observations to knots, in order to include the effect of the interaction between the covariates in the model. One of the advantages of this model is that the estimated smooth curve allows for the adjustment of a wide variety of nonlinear functions, allowing its application in a wide variety of real situations. The regularization of the model is obtained by penalizing the parameters via penalized maximum likelihood (PML) which has the advantage of producing a shrinkage of the coefficient estimates and reducing overfitting. These methods are equivalent to the optimization of the constrained maximum likelihood and also to Bayesian methods in which the coefficients have a priori normal distributions with zero mean. The deviance test was used to validate the fitted model. The results showed that the adjusted model was significantly better than the stationary model with a reliability of 99%.
Regarding the empirical analysis of UVB radiation on the metropolitan area of Mexico City, we characterized the distribution of maxima in the spatial and temporal plane. In the spatial plane, although the results show the existence of differentiated local patterns, the estimates of the location parameter of the GEV distribution showed that there is a plane that determines the trend in the entire study region, which evidences the existence of a positive linear correlation in the west direction of the study area. These results are consistent with the demographic characteristics of the area. In the temporal plane, we observe cyclical observations on the location parameter of the spatial distribution of maxima. Such oscillations are dominated by a negative linear trend with respect to time, which is consistent with the increase in population and its corresponding consequences on air pollution. Our findings also revealed the existence of areas with well-defined spatio-temporal patterns which should help administrative authorities to improve prevention policies and standards to mitigate the impact of UVB radiation maxima. Regarding the simulation results, these demonstrate that it is feasible to identify the nonlinear characteristics of the trends reliably under the parametric conditions used in the simulation, which were established with values of the GEV parameters similar to those found in real conditions. Particularly, the spatial function for the location parameter used in the simulation, which contains nonlinear features that we can expect to find in real data, was satisfactorily estimated. However, we conclude that optimal simulations related to the distribution of nonstationary GEV is an issue that requires further investigation.
Future work can first include to analyze new functions of distance between vectors, since it is natural to think that some variables may be more important as explanatory variables than others. Secondly, to examine the asymptotic properties of the estimators regarding to the number of knots. Finally, some other further studies could consider the analysis of the sensitivity of estimates on more complex nonlinear functions using simulations on the GEV distribution under different sample sizes.