Estimating Precipitation Extremes using Log-Histospline
Whitney K. Huang Douglas W. Nychka Hao Zhang










Estimating Precipitation Extremes using the Log-Histospline
    Whitney K. Huang[Statistical and Applied Mathematical Sciences Institute. E-mail: whuang@samsi.infowhuang@samsi.info and Pacific Climate Impacts Consortium, University of Victoria. E-mail: whitneyhuang@uvic.cawhitneyhuang@uvic.ca], Douglas W. Nychka[Department of Applied Mathematics and Statistics, Colorado School of Mines. E-mail: nychka@mines.conychka@mines.co], Hao Zhang[Department of Statistics, Purdue University. E-mail: zhanghao@purdue.eduzhanghao@purdue.edu]
    December 30, 2023
=========================================================================================================================================================================================================================================================================================================================================================================================================================================================================================




 
One of the commonly used approaches to modeling extremes is the peaks-over-threshold (POT) method. The POT method models exceedances over a threshold that is sufficiently high so that the exceedance has approximately a generalized Pareto distribution (GPD). This method requires the selection of a threshold that might affect the estimates. Here we propose an alternative method, the Log-Histospline (LHSpline), to explore modeling the tail behavior and the remainder of the density in one step using the full range of the data. LHSpline applies a smoothing spline model to a finely binned histogram of the log transformed data to estimate its log density. By construction, a LHSpline estimation is constrained to have  polynomial tail behavior, a feature commonly observed in daily rainfall observations. We illustrate the LHSpline method by analyzing precipitation data collected in Houston, Texas.















§ INTRODUCTION
 

Estimating extreme quantiles is crucial for risk management in a variety of applications. For example, an engineer would seek to estimate the magnitude of the flood event which is exceeded once in 100 years on average, the so-called 100-year return level, based on a few decades of time series <cit.>. A financial analyst needs to provide estimates of the Value at Risk (VaR) for a given portfolio, essentially the high quantiles of financial loss <cit.>. In climate change studies, as the research focus shifts from estimating the global mean state of climate variables to the understanding of regional and local climate extremes, there is a pressing need for better estimation of the magnitudes of extremes and their potential changes in a changing climate <cit.>.



The estimation of extreme quantiles poses a unique statistical challenge. Essentially, such an estimation pertains to the upper tail of a distribution where the available data of extreme values are usually sparse (e.g. see Fig. <ref>). As a result, the estimate will have a large variance that can increase rapidly as we move progressively to high quantiles. Furthermore, if the quantile being estimated is beyond the range of data (e.g., estimating the 100-year return level given the 50 years history of observation), the need to explicitly model the tail with some parametric form is unavoidable. Extreme value theory (EVT) provides a mathematical framework of performing inference for the upper tail of distributions. One widely used approach, based on the extremal types theorem <cit.>, is the so-called block maxima (BM) method where one fits a generalized extreme value (GEV) distribution  to block maxima given that the block size is large enough. The reader is referred to <cit.>, <cit.> and <cit.> for more details.




One drawback of the BM method is that data are not used efficiently, that is, the block size b needs to be sufficiently large so that the GEV distribution is approximately valid for a given time series (with sample size n), which makes the sample size of the block maxima (n/b) substantially smaller than the original sample size. Moreover, the top few largest order statistics within a given block should, in principle, inform us about the behavior of extreme events <cit.>. The peak-over-threshold (POT) method is based on the Pickands–Balkema–de Haan theory <cit.>. For a random variable Y with cumulative distribution function F(·), it states that the distribution of the exceedance over the threshold (i.e., the distribution of Y-u given Y>u) converges to a generalized Pareto distribution (GPD) as the threshold u tends to the upper limit point y_F = sup{ y: F(y) < 1}. With an appropriately chosen u, the POT method makes use of the exceedances and the censored data of those that do not exceed the threshold <cit.>. One advantage of the POT method over the BM method is that it typically makes use of the available data more efficiently in estimating extreme events.



However, to apply the POT method, a threshold u has to be chosen and the estimates may be sensitive to the chosen threshold <cit.>. The threshold selection involves a bias-variance trade-off: if the chosen threshold is too high, the estimates exhibit large variability; if the chosen threshold is too low, the asymptotic justification of the GPD approximation to the tail density is less effective, which leads to bias. There are several graphical tools that aim to help with the threshold selection, for example, the mean residual life plot and the parameter stability plot <cit.>. However, the use of these graphical tools does not always lead to a clear choice of the threshold. In general, automated threshold selection is a difficult problem <cit.>. In addition, the POT method does not assume anything about the distribution below the chosen threshold, which can be of interest in some applications (e.g., stochastic weather generators <cit.>), although some efforts have been made. We will review some attempts in the next paragraph.




Recently, there are some attempts to model the distribution of a random variable while retaining a GPD tail behavior <cit.>. These methods can be broadly divided into two categories: the extreme value mixture models  and the extended generalized Pareto distribution (EGPD) method <cit.>. The basic idea of the extreme value mixture is to model a distribution as a mixture of a bulk distribution and a GPD distribution above a threshold. The threshold is then treated as an additional parameter to be estimated from the data. However, the specification of the bulk distribution can have non-negligible impacts on the estimates of GPD parameters, which can lead to substantial biases in the tail estimates. Finally, the estimation uncertainty for the threshold can be quite large <cit.> and hence it is not clear whether this parameter can be identifiable. The extended Pareto method proposed by <cit.> bypasses this issue of the mixture modeling by proposing several classes of parametric models with GPD limiting behavior for the upper tail. <cit.> extended the scope of this approach by allowing the classes of parametric models with GPD limiting behavior for both lower and upper tails in an application to rainfall modeling.



This work presents a new statistical method, called the Log-Histospline (LHSpline), for estimating probabilities associated with extreme values. Similar to <cit.>, this method applies  extreme value theory to the upper tail distribution. However, it does not impose a parametric form in the bulk of the distribution where the density can be determined from the data. This work is motivated by some applications in climate studies, one of which involves precipitation data. We would like to (i) provide a flexible model to the full range of the non-zero rainfall distribution, and (ii) reliably estimate extreme events (e.g. 100-year daily rainfall amount).  Specifically, our approach is to first log transform the nonzero daily rainfall observations and then apply a generalized cubic smoothing spline on a finely binned histogram to estimate the log-density. The purpose of the data transformation step is similar to that of a variable bandwidth in  density–estimation literature <cit.>. Applying the spline smoothing log-density estimation <cit.> to the transformed variable will ensure the algebraic (e.g., Pareto) upper tail behavior and gamma-like lower (near 0) tail behavior, commonly observed in daily rainfall data.



The rest of this paper is structured as follows: In Section <ref>, the LHSpline method is introduced and computation for inferences is described. A simulation study is presented in Section <ref>. An application to daily precipitation collected in Houston Hobby Airport, Texas is illustrated in Section <ref>. Some discussion of future research directions is provided in the last section.   




§ THE LOG-HISTOSPLINE
 



 §.§ The Model of Log Density: Natural Cubic Spline

Let Y be a continuous random variable with a probability density function (pdf) f(y). Throughout this work we assume that Y ∈ (0, ∞] is heavy-tailed, i.e.,

    f(y) ∼ y^-(α+1) as  y →∞, α >0.
  

To represent the heavy tail of the distribution of Y,  we take the logarithmic transformation X = log (Y), via its log-density, and assume the pdf of X takes the following form:

    e^g(x),    x ∈ (-∞, ∞]
  
Specifically, 
the density f(·) and log-density g(·) are related in the following way:

    f(y) =  y^-1e^g(log(y)),       y > 0.
 

Modeling the log-density g has an advantage that it conveniently enforces the positivity constraint on f
<cit.>. Here we would like to model g nonparametrically to avoid a strong and likely misspecified parametric structure <cit.>. However, it is well-known that the usual nonparametric density estimation procedure will often introduce spurious features in the tail density due to limited data in the tail <cit.>. This issue is amplified when dealing with heavy-tailed distributions that are typically observed in precipitation data. On the other hand, since the focus is on estimating extreme quantiles, it is critical that the model can capture both qualitative (polynomial decay) and quantitative (the tail index α) behaviors of the tail density. Here we assume that g belongs to the family of natural cubic splines (see Definition <ref>) to accommodate both a flexible bulk distribution and the algebraic tail. The LHSpline combines the ideas of spline smoothing <cit.> and Histospline <cit.> to derive an estimator to achieve these goals.
To facilitate the description of our method, we review the following definition <cit.>
 
A natural cubic spline g(x), x ∈ [a,b] is defined with a set of points {ζ_i}_i=1^N (knots) such that -∞≤ a ≤ζ_1≤ζ_2≤⋯≤ζ_N < b ≤∞ with the following properties:


  * g is linear, x ∈ [a, ζ_1],    x ∈ [ζ_N, b]

  * g is a cubic polynomial, x ∈ [ζ_i, ζ_i+1],     i = 1, ⋯, N -1

  * g ∈𝒞^2,     x∈ (a, b)
   
where  𝒞^k is the class of functions with k continuous derivatives 


The cubic spline smoother is widely used in nonparametric density estimation due to its flexibility  and hence provides a flexible model for the bulk of the distribution <cit.>. 
In our setting,  when y ≥exp(ζ_N), the tail density of Y is  

    f(y)     = y^-1e^g(log(y)) = y^-1e^-αlog(y) + β
        = y^-1 e^-αlog(y)e^β = Cy^-1e^-αlog(y)
        = Cy^-1y^-α = Cy^-(α + 1)4

where -α is the slope and β is the intercept of the log-density g at ζ_N. Therefore, by the linear boundary conditions on a cubic spline we automatically obtain the algebraic right tail behavior. 



 §.§ The Estimation Procedure
 

Given an i.i.d. sample Y_1, Y_2, …, Y_n of Y, the estimation of g involves two steps: data binning to create a histogram object, and a generalized smoothing spline of the histogram. In what follows, it is useful (see Fig. <ref>) to illustrate these steps using a synthetic data set simulated from a EGPD . Some considerations when applying LHSpline to precipitation data will be discussed in Section <ref>.



Data binning: First, bin the transformed data X_i=log(Y_i), i=1, …, n by choosing N+1 equally spaced break points  {b_j } to construct the corresponding histogram object (i.e. the histogram counts {Z_j} = #{ X_i: X_i∈ [b_j, b_j+1]}. 
We set the knots associated with the spline to the bin midpoints:
ζ_j = (b_j + b_j+1)/2. Several remarks on data binning should be made here:


  * Equal-sized bins in the log scale implies a variable bin size with the bin size becoming increasingly larger in the upper tail.

  * The number of bins should be “sufficiently” large for the Poisson assumption made in the next step justifiable.

  * The choice of the first and especially the last break points (i.e. b_1 and b_N+1) will have non-negligible impacts on tail estimation in our framework. We will choose them somehow smaller (larger) than the minimum X_(1) (maximum X_(n)) of the log-transformed data, which is different than what is typically done in constructing a histogram. We will defer the discussion on this point to Sec <ref>.   

  * It is assumed  that the bin size is fine enough so that  ∫_b_i^b_i+1 e^g(x) dx is well approximated by 
e^g(ζ_i)(b_i+1 - b_i).


Smoothing the histogram: We adapt a penalized approach <cit.> to obtain a functional estimate of g as follows. First, we assume the bin counts {Z_j}_j=1^N each follow a Poisson distribution[{Z_j}_j=1^N has a multinomial distribution with parameters n and π = (π_1, ⋯, π_N) where π_j is the probability that a random variable X falls into the j_th bin. Here we consider Z_jind∼Pois(γπ_j),   j = 1,⋯, N. The MLE of γ under the Poisson model follows γ̂ = ∑_j=1^NZ_j = n and π̂ is equal to the MLE of π under the multinomial model.] <cit.> with log intensity g̃_j,   j = 1, ⋯, N, given that N is large enough so that the bin size is fine enough. Furthermore, we assume g̃_j≈g̃(ζ_j), i.e., we assume the Poisson log intensity at each bin represents the log intensity function (unnormalized log density function) evaluated at its midpoint. Hence we perform a penalized Poisson regression to the data pair {Z_j, ζ_j}_j=1^N using a log link function and with a penalty term that penalizes the “roughness” of g̃. The estimate g̃, the unnormalized version of g, can be obtained by solving the following minimization problem:


    -∑_j=1^N{g̃(x_j)z_j-e^g̃(x_j) -log(z_j)!
    } 
    + λ(∫ _x ∈ℝ(g̃^''(x))^2 dx)

where λ is the smoothing parameter. Note that the first term in the objective function is the negative log likelihood for the histogram counts under a Poisson approximation to the distribution, and the second term, the squared integral of the second derivative of g̃ multiplied by the smoothing parameter, is the penalty function. The solution to this optimization problem exists, is unique, and has the form of a natural cubic spline for any λ <cit.>. The selection of the smoothing parameter λ plays a role of balancing the data fidelity of the Poisson regression fit to the histogram counts, as presented by the negative log-likelihood, and the “smoothness” of g̃, and is chosen by approximate cross validation . Note that the smoothness penalty favors linear functions, that correspond to the algebraic tail behavior we want in the untransformed distribution. Finally we renormalize the g̃ to make the integral of ∫_ℝ e^g equal to one. Numerical integration is used to approximate the integral within the data range and an analytic form is used for the density beyond the bin limits.






 §.§ Bias Correction
 



  §.§.§ Extending Data Range
 

As mentioned in the remark (3) in the data binning step, the choice of b_1 and b_N+1 plays an important role in our penalized approach. Under the usual histogram construction (i.e. b_1 = X_(1), b_N+1 = X_(n), see Fig <ref>) and a large number of bins, as in our setting, one will very likely observe “bumps” (i.e., isolated non-zero counts) in the boundaries of the histogram. As a result, the cross validation will maintain the smoothness and tend to overestimate the slopes near the boundary knots. Our solution to reduce this bias is to extend the range of the histogram beyond the sample maximum/minimum so that some zero counts beyond the boundaries will be included. In fact, one can think of the estimation procedure of the LHSpline as a discrete approximation to estimating the intensity (density) function of an  (normalized) inhomogeneous Poisson process <cit.>. In this regard, one should take into account the support (observation window in the point process context) as part of the data. 



  §.§.§ Further Corrections: Bootstrap and Smoothing Parameter Adjustment
 

After applying the aforementioned boundary correction by extending the observation window, tail bias still exists.  Here we propose two approaches to correct the bias. The first approach is to estimate the bias via a “parametric”  bootstrap <cit.> on the LHSpline estimate. Note that given the approximate Poisson model for the histogram bin counts it is straightforward to generate bootstrap samples for this computation. The bias is estimated as the difference between the (point-wise) mean of the bootstrap estimates g̅^̅*̅ and the “true” g, in this case, ĝ estimated from the histogram.  
We then subtract this bias term from the LHSpline estimate.  The second approach is to adjust the smoothing parameter estimate λ̂_CV directly. We found that the smoothing parameter obtained from CV may not give the best result as here the primary objective is to estimate the tail density. It was found that, in our LHSpline, reducing λ,   0.05 ×λ̂_CV,  gives a better tail estimation while still maintaining a good estimation performance in the bulk distribution (see Sec <ref>).    



 §.§ Quantifying Estimation Uncertainty: The “Bayesian” Approach
 

Here we give a brief account of the Bayesian interpretation of smoothing splines that will be used for quantifying estimation uncertainty for the LHSpline. Readers are referred to  for more details.


From a Bayesian perspective, the penalty term in Equation (<ref>) effectively introduces a Gaussian process prior on g, the log density. Hence the LHSpline estimation can be thought as a procedure of finding the posterior mode of g with a Poisson likelihood for the histogram and a zero mean Gaussian process with a generalized covariance (the reproducing kernel), K(·,·), as the prior. The LHSpline estimate of a given histogram can be found by using the Fisher scoring algorithm  where the minimization of the non-quadratic negative penalized log-likelihood is approximated by the method of iteratively reweighted least squares (IRLS).


The goal here is to sample from an approximate posterior distribution [g|u, λ] where g = {g _k = g(ζ_k)}_k=1^N, u = {u_k}_k=1^N are the “pseudo” observations in IRLS, λ is the smoothing parameter, and the likelihood has been approximated with a multivariate Gaussian distribution. Based on the linear approximation around the posterior model we have

    [g|u, λ] ∼MVN(ĝ, (W + Γ)^-1)

where W is a diagonal matrix representing the approximate precision matrix of the pseudo observations and Γ is the prior precision for g based on the Bayesian interpretation of a spline. Several assumptions are being made here: 1) the smoothing parameter λ is known; 2) the linear problem is assumed to be multivariate normal and finally 3) a proper prior, N(0, ϵ I), for the linear trend (i.e., the null space of a natural cubic spline) has been taken in the limit to be improper (i.e., ϵ→∞).


The sampling approach is similar to a classical bootstrap <cit.> with the following steps: 


  * generate g^*∼MVN(0, Γ^-1)

  * generate pseudo observations u^*∼MVN(g^*, W^-1)

  * compute the estimate based on u^*: ĝ^* = (W+Γ)^-1Wu^*

  * compute error: u = g^* - ĝ^*

  * approximate draw from the posterior, ĝ + u, where ĝ is the LHSpline estimate (the posterior mode/mean under the aforementioned Gaussian assumption) evaluated at knots {ζ_k}_k=1^N. 





§ SIMULATION STUDY
 

The purposes of this simulation study are threefold: (i) to demonstrate how we implement the LHSpline, (ii) to compare the LHSpline method with the EVT-based methods (i.e. POT and BM) in estimating extreme quantiles, and (iii) to compare with a gamma distribution fit to the bulk distribution.




 §.§ Setup
 
We conduct a Monte Carlo study by simulating the data from the model (i) of the EGPD . The basic idea of EGPD is to modify the GPD random number generator H_σ, ξ^-1(U), where H^-1_σ, ξ denotes the inverse GPD cdf with scale (σ) and shape (ξ) parameters, by replacing U, the standard uniform random variable, with V = G^-1(U) where G is a continuous cdf on [0,1]. The model (i) takes G(v) = v^κ, κ>0 which contains the GPD as a special case (i.e., κ = 1).


We choose the parameter values to be κ = 0.8, σ = 8.5, ξ =0.2 and the sample size n = 18,250, which corresponds to a time series with 50 years of complete daily data (ignore the leap years). The EGPD parameters are chosen to reflect typical distributional behavior of daily precipitation. We repeat the Monte Carlo experiment 100 times and evaluate the estimation performance for 25-, 50-, and 100- year return levels.




 §.§ LHSpline Illustration
 
We use a simulated data series (see Fig. <ref>) to illustrate the LHSpline method as described in Sec. <ref>. To illustrate the “boundary effect” we first set equidistant breakpoints so that the knots (the midpoints of the histogram) span the whole range of the data (i.e. ζ_1 = x_(1), ζ_N=x_(n) where ζ_1 and ζ_N are the first and last knot points). We choose N = 150 meaning that we construct a histogram with a rather large number of bins (much larger than what the default in  in  would suggest). We then apply the generalized smoothing spline procedure by solving equation (<ref>) via approximate cross-validation to obtain an estimate. Figure  <ref> shows the LHSpline estimate along with a GPD fit (with threshold u chosen as the .95 empirical quantile) and the true density. Upon visual inspection one may conclude that LHSpline performs well, or at least as good as the POT method, and both estimates are fairly close to the true density. 

 

However, a more careful examination for the log-log plot (log-density against log(y), see the left panel of Fig <ref>) and the return levels estimation reveals that LHSpline with ζ_1 = x_(1) and ζ_N=x_(n), in general, overestimates the extreme quantiles (see Fig. <ref>, right panel).

 

The issue of overestimation of LHSpline is partly due to the “boundary effect” when fitting penalized Poisson regression to a histogram. Specifically, the histogram counts in Fig. <ref> are {Z_1 = 1, Z_2 = Z_3 = ⋯ = Z_26 = 0, ⋯,Z_148 = Z_149 = 0, Z_150 = 1}. The histogram counts in the first and the last bins are non-zero (typically 1) by construction but the nearby bins are likely to be zero. Therefore, these “bumps” at both ends force the smoothing spline to overestimate the slopes at the boundaries (sometimes it will give nonsensical results, e.g., negative slope at the left boundary or/and positive slope at the right boundary) and hence the extreme quantiles.
 

To alleviate this boundary effect, we extend the range and number of bins in the histogram beyond the range of the data to include extra bins that have zero counts and refit these augmented data pairs (the original counts and those extra zero counts) to penalized Poisson regression. We observe that this strategy does remove some of the positive bias in return levels estimation (see Fig. <ref>) but there is still some positive bias remaining. We then apply the bias correction via bootstrap with 200 bootstrap samples and the aforementioned smoothing parameter adjustment. Both approaches further reduce the bias (see Fig <ref>). 

 





 §.§ Estimator Performance
 

In this subsection we first assess the performance of estimating high quantiles using an LHSpline and two commonly used EVT-based methods, namely, the block maxima (BM) method by fitting a GEV to “annual maxima”, and the peaks-over-threshold (POT) method by fitting a GPD to excesses over a high threshold (.95 empirical quantile in our study). In order to put this comparison on an equal footing as much as possible, we put a positive Gamma prior (with mean equal to 0.2, the true value) on ξ
when fitting the GEV and GPD models.

We use the  function in the 
R package <cit.> with the  Generalized maximum likelihood estimation (GMLE) method <cit.> to estimate the GEV and GPD parameters.
We compare the GMLE plug-in estimates for 25-, 50-, and 100-year return levels with that of the LHSpline estimates (with and without bias correction). To get a better sense of how well each method performs, we also include the “oracle estimates” where we fit the true model (EGPD) using the maximum likelihood method to obtain the corresponding plug-in estimates.
 

Figure <ref> shows that, perhaps not surprisingly, the variability of the return level estimates obtained by fitting a GEV distribution to block maxima is generally larger than that of the estimates obtained by fitting a GPD distribution to threshold exceedances under the i.i.d. setting here. Also, it is clear that the LHSpline without bias correction can not only lead to serious bias but also inflate the estimation variance. After applying two different bias corrections mentioned in Sec <ref>, the estimate becomes closer to being unbiased (especially the bias correction with adjustment of the smoothing parameter) with a reduction in the estimator variability. Also the estimator performance is comparable with that of the GPD.



We then assess the estimation uncertainty using the conditional simulation approach described in Sec <ref> and compare with the GPD interval estimation obtained by the delta and profile likelihood methods. To simplify this presentation, we only show the result of the 90% confidence interval for the 50-year return level (see Fig <ref>)




We summarize in Table <ref> the empirical coverage probability (ECP) and the width of CI for the  delta and profile likelihood methods when fitting a GPD to excesses over the .95 quantile and the conditional simulation (cond-sim) when fitting LHSpline with λ = 0.05 ×λ̂_CV for 25-, 50-, and 100-year return levels, respectively. 




We also assess the estimation performance of LHSpline in the bulk distribution. Fig <ref> shows that the LHSpline performs much better than the gamma fit, a widely used distribution for modeling precipitation <cit.>, and performs nearly as well as the “oracle” (EGPD) approach.






§ APPLICATIONS
 



 §.§ Motivation: Hurricane Harvey Extreme Rainfall
 

Hurricane Harvey brought unprecedented amounts of rainfall to the Houston metropolitan area between the 25th and the 31st of August 2017, resulting in catastrophic damages to personal property and infrastructure. In the wake of such an extreme event, there is  interest in understanding its rarity and how human-induced
climate change might alter the chances of observing such an event <cit.>. In this section, we apply the LHSpline to the daily precipitation measurements for Houston Hobby airport (see Fig. <ref> for the histogram and Fig. <ref> for the time series) from the Global Historical  Climatology Network (GHCN) <cit.>.



 §.§ Houston Hobby Rainfall Data Analysis
 

We fit a LHSpline with 150 equally spaced bins to the full range of non-zero precipitation for Hobby Airport (∼ 27.6% of all the observations) prior to this event (from Jan. 1949 to Dec. 2016). To facilitate a comparison with the POT method we also fit a GPD to excesses over a high threshold (chosen as u_0 = 43  mm ∼ .93 empirical quantile of nonzero daily precipitation, see diagnostic plot in Fig. <ref>).
  







An important practical issue when applying the LHSpline is that the discretization of rainfall measurement (1/100^th of an inch) introduces an artifact in the finely binned histogram in the log scale (see Fig. <ref>, left panel) and hence the cross validation choice for λ is affected by this discretization. In practice, all the precipitation measurements have precision limits and thus it is important to take this fact into account especially when the measurements are small (i.e., near zero). Here we treat the data being left-censored by truncating the values below two different values (exp(1) and exp(2)) to alleviate this “discreteness” effect.



As has been demonstrated in Fig. <ref> and Fig. <ref> in Sec. <ref>, a good visual agreement in log-density (density) does not necessarily imply that it actually provides a good return level estimate. Fig. <ref> shows the estimated log-log plot and return levels ranging from 2 years to 100 years (under the assumption of temporal stationarity). The LHSpline gives somewhat lower estimates than that of the GPD estimates. A quantile-quantile plot in Fig. <ref> (left panel) indicates that there might be an issue of overestimating extreme high quantiles in the GPD fit. That is, the GPD estimate of the 68-year rainfall is 395.47 (90% CI (300.88, 569.67)) mm whereas the (λ-adjusted) LHSpline estimates (with lower bound 0.254, exp(1), and exp(2)) are 371.04 (302.39, 567.62) mm, 310.64 (261.62, 491.80) mm, and 287.06 (246.86, 457.05) mm, respectively, which are somewhat closer to the maximum value (252.73 mm) during the 1949 ∼ 2016 period. However, one should be aware that the sample maximum can be quite variable and hence might not reflect the true magnitude of the “68-year rainfall”. Fig. <ref> (right panel) also indicates that the LHSpline might provide a better fit than the EGPD to the Hobby Airport data.    



 

    < g r a p h i c s >

The quantile-quantile plot of the Hobby daily precipitation data. Left: A comparison between LHSplines and GPD for the upper 5% (≥ 50  mm) of rainy data). Right: A comparison between LHSplines and EGPDs for the “full range” (left truncated at exp(1).



Lastly, we would like to investigate the question: “How unusual was the event of 306.58 mm (26th August, 2017) in daily total precipitation at Hobby Airport?” To simplify the matter we make a stationarity assumption, that is, the distribution of daily precipitation does not change during the time period 1949 ∼ 2017. The estimates for the return period of this event are given in the following table:



The much shorter return period estimate obtained from the POT method is due to the overestimation of the upper tail (see Fig. <ref>) whereas the somewhat longer return periods estimated by LHSplines might be more aligned with what people would expect. Although again one should be aware that there exists, among other issues, a large estimation uncertainty with respect to large return periods (see Table. <ref>).  





§ DISCUSSION
 

This work presents a new statistical method, the LHspline, for estimating extreme quantiles of heavy-tailed distributions. In contrast with some widely used EVT based methods that require extracting “extreme” observations to fit a corresponding asymptotically justified distribution, the LHSpline makes use of the full range of the observations for the fitting. By combining data transform and spline smoothing, the LHSpline estimation effectively achieved the desirable tail structure that is consistent with EVT and a flexible bulk distribution in the context of rainfall modeling. We demonstrate through simulation that this method
performs comparable to the POT method for return level estimation
with an additional benefit that it jointly models the bulk and the tail of a distribution.


However, by construction, the LHSpline method is only applicable for heavy-tailed distributions, which excludes many important environmental processes, for example, air surface temperature which may have a bounded tail <cit.>. In terms of implementation, the LHSpline requires some additional tuning, for example, the number of bins for the histogram should grow with the sample size. Limited experiments (not reported here) suggest the estimate is not sensitive to the number of bins once the binning is chosen fine enough. Another tuning issue is to decide how far one should extend beyond the range of the observations and how to adjust the smoothing parameter λ to remove the boundary effect. Here we suggest extending the bin range and reducing the smoothing parameter. This choice is in place of picking the threshold in the POT approach and we believe it is less sensitive in terms of density estimation.


The theoretical properties of our method are largely unexplored. Much of the theoretical results developed in nonparametric density estimation concern the performance in terms of global indices such as 𝔼[∫_x ∈ℝ |f(x) - f̂(x)|  dx] or 𝔼[sup_x ∈ℝ |f(x) - f̂(x)|] and are largely confined in a bounded region (i.e. x ∈ [a,b], a and b finite). It is not clear how these results can inform us about the estimation performance for extreme quantiles on a potentially unbounded domain.




Applying LHSpline to many observational records or the grid cells of high resolution and/or ensemble climate model experiments <cit.> will result in summaries of the distribution that are well suited for further data mining and analysis. The log-density form is particularly convenient for dimension reduction because linear (additive) projection methods such as principle component analysis make sense in the log space.








In general we believe the LHSpline will be a useful tool as climate informatics tackles complex problems of quantifying climate extremes. 



§ ACKNOWLEDGEMENTS
 
This work was conducted as part of the https://www.statmos.washington.eduResearch Network for Statistical Methods for Atmospheric and Oceanic Sciences (STATMOS), supported by the National Science Foundation (NSF) Grants #s 1106862, 1106974, and 1107046. This work was also supported by the NSF Grant # 1638521 to the Statistical and Applied Mathematical Sciences Institute(SAMSI), the NSF Grant # 1406536, and by the National Center for Atmospheric Research (NCAR). NCAR is sponsored by the NSF and managed by the University Corporation for Atmospheric Research. The authors would like to thank the anonymous reviewers for their valuable comments and suggestions to improve the quality of the paper.





