Modeling Precipitation Extremes using Log-Histospline

One of the commonly used approaches to modeling univariate extremes is the peaks-over-threshold (POT) method. The POT method models exceedances over a (sufficiently high/low) threshold as 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 geophysical observations. We illustrate LHSpline method by analyzing precipitation data collected in Houston, Texas.



page 1

page 2

page 3

page 4


Functional Peaks-over-threshold Analysis

Peaks-over-threshold analysis using the generalized Pareto distribution ...

Trimming and threshold selection in extremes

We consider removing lower order statistics from the classical Hill esti...

On a minimum distance procedure for threshold selection in tail analysis

Power-law distributions have been widely observed in different areas of ...

Modelling extreme claims via composite models and threshold selection methods

The existence of large and extreme claims of a non-life insurance portfo...

The application of accumulation tests in Peaks-Over-Threshold modeling with Norwegian Fire insurance Data

Modeling excess remains to be an important topic in insurance data model...

Estimation of Extreme Survival Probabilities with Cox Model

We propose an extension of the regular Cox's proportional hazards model ...

Modeling Variables with a Detection Limit using a Truncated Normal Distribution with Censoring

When data are collected subject to a detection limit, observations below...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 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

(Katz et al., 2002). A financial analyst needs to provide estimates of the Value at Risk (VaR) for a given portfolio, essentially the high quantiles of financial loss (Embrechts et al., 1997; Tsay, 2005). 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 (Zwiers and Kharin, 1998; Easterling et al., 2000; Alexander et al., 2006; Tebaldi et al., 2006; Cooley et al., 2007; Kharin et al., 2007; AghaKouchak et al., 2012; Shaby and Reich, 2012; Huang et al., 2016).

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. 1

). 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 (Fisher and Tippett, 1928; Gnedenko, 1943), 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 (Jenkinson, 1955), Gumbel (1958) and Coles (2001) for more details.

Figure 1:

Histogram of daily rainfall amount (mm) at the Hobby Airport in Houston, Texas, from 1949 to 2017. The vertical ticks at the x-axis are the values of the individual data points and the black dashed line is a kernel density estimate. Due to its tail heaviness, the largest values are substantially larger than the bulk of the distribution. Three out of the five largest precipitation measurements (indicated by red arrows) were observed during Hurricane Harvey, August 25 - 31, 2017

One drawback of the BM method is that data are not used efficiently, that is, the block size needs to be sufficiently large so that the GEV distribution is approximately valid for a given time series (with sample size ), which makes the sample size of the block maxima 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 (Weissman, 1978; Smith, 1986). The peak-over-threshold (POT) method is based on the Pickands–Balkema–de Haan theory (Pickands, 1975; Balkema and De Haan, 1974)

. For a random variable

with cumulative distribution function

, it states that the distribution of the exceedance over the threshold (i.e., the distribution of given ) converges to a generalized Pareto distribution (GPD) as the threshold tends to the upper limit point . With an appropriately chosen , the POT method makes use of the exceedances and the censored data of those that do not exceed the threshold (Smith, 1984; Davison, 1984; Davison and Smith, 1990). 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 has to be chosen and the estimates may be sensitive to the chosen threshold (Scarrott and MacDonald, 2012; Wadsworth and Tawn, 2012). 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 (e.g. Coles, 2001, p.80 and p.85). 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 (Dupuis, 1999; Wadsworth and Tawn, 2012). 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 (Kleiber et al., 2012; Li et al., 2012)), 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 (Frigessi et al., 2002; Tancredi et al., 2006; Carreau and Bengio, 2009; Papastathopoulos and Tawn, 2013; Naveau et al., 2016). These methods can be broadly divided into two categories: the extreme value mixture models see Scarrott and MacDonald, 2012; Dey and Yan, 2016 for review and the extended generalized Pareto distribution (EGPD) method (Papastathopoulos and Tawn, 2013; Naveau et al., 2016). 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 (Frigessi et al., 2002) and hence it is not clear whether this parameter can be identifiable. The extended Pareto method proposed by Papastathopoulos and Tawn (2013)

bypasses this issue of the mixture modeling by proposing several classes of parametric models with GPD limiting behavior for the upper tail.

Naveau et al. (2016) 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

Naveau et al. (2016), 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 (Wand et al., 1991). Applying the spline smoothing log-density estimation (Silverman, 1982; Gu and Qiu, 1993) 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 2, the LHSpline method is introduced and computation for inferences is described. A simulation study is presented in Section 3. An application to daily precipitation collected in Houston Hobby Airport, Texas is illustrated in Section 4. Some discussion of future research directions is provided in the last section.

2 The Log-Histospline

2.1 The Model of Log Density: Natural Cubic Spline


be a continuous random variable with a probability density function (pdf)

. Throughout this work we assume that is heavy-tailed, i.e.,


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


Specifically, the density and log-density are related in the following way:


Modeling the log-density has an advantage that it conveniently enforces the positivity constraint on (Leonard, 1978; Silverman, 1982; Kooperberg and Stone, 1991; Eilers and Marx, 1996). Here we would like to model nonparametrically to avoid a strong and likely misspecified parametric structure (e.g. Silverman, 1986). 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 (e.g. Kooperberg and Stone, 1991). 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 belongs to the family of natural cubic splines (see Definition 1) to accommodate both a flexible bulk distribution and the algebraic tail. The LHSpline combines the ideas of spline smoothing (e.g. Wahba, 1990; Gu, 2013) and Histospline (Boneva et al., 1971) to derive an estimator to achieve these goals. To facilitate the description of our method, we review the following definition (Wahba, 1990)

Definition 1 (Natural cubic spline)

A natural cubic spline , is defined with a set of points (knots) such that with the following properties:

  1. is linear,

  2. is a cubic polynomial,

where is the class of functions with 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 (e.g. Silverman, 1986; Gu, 2013). In our setting, when , the tail density of is


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

2.2 The Estimation Procedure

Given an i.i.d. sample of , the estimation of 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. 2) to illustrate these steps using a synthetic data set simulated from a EGPD Naveau et al., 2016, p. 2757. Some considerations when applying LHSpline to precipitation data will be discussed in Section 4.

Data binning: First, bin the transformed data by choosing equally spaced break points to construct the corresponding histogram object (i.e. the histogram counts . We set the knots associated with the spline to the bin midpoints: . Several remarks on data binning should be made here:

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

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

  3. The choice of the first and especially the last break points (i.e.  and ) will have non-negligible impacts on tail estimation in our framework. We will choose them somehow smaller (larger) than the minimum (maximum ) 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 2.3.

  4. It is assumed that the bin size is fine enough so that is well approximated by

Smoothing the histogram: We adapt a penalized approach (Good and Gaskins, 1971; Tapia, 1978; Green and Silverman, 1994) to obtain a functional estimate of as follows. First, we assume the bin counts

each follow a Poisson distribution

444{ has a multinomial distribution with parameters and where is the probability that a random variable falls into the j bin. Here we consider . The MLE of under the Poisson model follows and is equal to the MLE of under the multinomial model. (Lindsey, 1974a, b; Eilers and Marx, 1996; Efron and Tibshirani, 1996) with log intensity , given that is large enough so that the bin size is fine enough. Furthermore, we assume , 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 using a log link function and with a penalty term that penalizes the “roughness” of . The estimate , the unnormalized version of , can be obtained by solving the following minimization problem:


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 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 (Green and Silverman, 1994; Gu, 2013). 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 , and is chosen by approximate cross validation CV, see Sullivan, 1988 for more details. 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 to make the integral of 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.

Figure 2: An illustration of applying LHSpline to a simulated data set. Here we bin the log transformed data into equally-spaced intervals to obtain a histogram. Red asterisks denote the non-zero counts and gray asterisks denote the zero counts. The LHSpline estimate (blue dashed line) is obtained by fitting a penalized Poisson regression to the histogram counts (both red and gray asterisks); the true density (the density of where ) is plotted as black solid line.

2.3 Bias Correction

2.3.1 Extending Data Range

As mentioned in the remark (3) in the data binning step, the choice of and plays an important role in our penalized approach. Under the usual histogram construction (i.e. , , see Fig 2) 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 (Brown et al., 2004). In this regard, one should take into account the support (observation window in the point process context) as part of the data.

2.3.2 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 (Efron, 1979) 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 and the “true” , 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 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 ,  , gives a better tail estimation while still maintaining a good estimation performance in the bulk distribution (see Sec 3.3).

2.4 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 Wahba, 1990, Chapter 5, Green and Silverman, 1994, Chapter 3.8, and Gu, 2013, Chapter 2.4 for more details.

From a Bayesian perspective, the penalty term in Equation (4) effectively introduces a Gaussian process prior on , the log density. Hence the LHSpline estimation can be thought as a procedure of finding the posterior mode of with a Poisson likelihood for the histogram and a zero mean Gaussian process with a generalized covariance (the reproducing kernel), , as the prior. The LHSpline estimate of a given histogram can be found by using the Fisher scoring algorithm Green and Silverman, 1994, p. 100 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 where , 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


where is a diagonal matrix representing the approximate precision matrix of the pseudo observations and is the prior precision for 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, , 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 (Efron, 1979) with the following steps:

  1. generate

  2. generate pseudo observations

  3. compute the estimate based on :

  4. compute error:

  5. approximate draw from the posterior, , where is the LHSpline estimate (the posterior mode/mean under the aforementioned Gaussian assumption) evaluated at knots .

3 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.

3.1 Setup

We conduct a Monte Carlo study by simulating the data from the model (i) of the EGPD see Naveau et al., 2016, p. 2757. The basic idea of EGPD is to modify the GPD random number generator , where denotes the inverse GPD cdf with scale and shape parameters, by replacing , the standard uniform random variable, with where is a continuous cdf on . The model (i) takes which contains the GPD as a special case (i.e., ).

We choose the parameter values to be and the sample size , 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 , , and year return levels.

3.2 LHSpline Illustration

We use a simulated data series (see Fig. 2) to illustrate the LHSpline method as described in Sec. 2. 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. , where and are the first and last knot points). We choose meaning that we construct a histogram with a rather large number of bins (much larger than what the default in hist in R would suggest). We then apply the generalized smoothing spline procedure by solving equation (4) via approximate cross-validation to obtain an estimate. Figure  3 shows the LHSpline estimate along with a GPD fit (with threshold chosen as the 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.

Figure 3: The density estimates for log scale (left) and original scale (right). In each panel, the blue dashed line is the LHSpline estimate (without boundary correction), red dot-dash POT estimate with as the .95 empirical quantile, and black truth density function.

However, a more careful examination for the log-log plot (log-density against , see the left panel of Fig 4) and the return levels estimation reveals that LHSpline with and , in general, overestimates the extreme quantiles (see Fig. 4, right panel).

Figure 4: Left: Estimated log density against . Right: Estimated return levels with return period ranging from 2 years to 100 years. The light blue dashed (red dot-dash) curves are 20 samples from the 100 Monte Carlo experiments for GPD fits with threshold quantile and LHSpline fits without boundary correction.

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. 2 are . 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. 5) 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 5).

Figure 5: The estimated log density and return levels with (green dashed line, LHS-er) and without (blue dashed line, LHS) boundary correction. Brown (purple) dashed lines is the LHSpline estimate with bootstrap ( adjustment) bias correction. The GPD estimate is shown in the red dot-dash line. The vertical red line is the chosen threshold for fitting the GPD.

3.3 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 ( 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 fevd function in the extRemes R package (Gilleland and Katz, 2016) with the Generalized maximum likelihood estimation (GMLE) method (Martins and Stedinger, 2000, 2001) 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 6 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 2.3, 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.

Figure 6: Boxplots of estimated Left: 25-, Middle: 50-, and Right: 100-year return levels using (from left to right in each panel) LHSpline1: LHSpline with boundary correction (extend data range by a factor of 1.5), LHSpline2: LHSpline with boundary correction and bootstrap bias correction, LHSpline3: LHSpline with boundary correction and smoothing parameter adjustment , GEV: block maxima method, GPD: Peaks-over-threshold method, and EGPD the oracle (extended GP). Boxplots are based on 100 independent replicates, and true values are represented by horizontal red lines.

We then assess the estimation uncertainty using the conditional simulation approach described in Sec 2.4

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 


Figure 7: The interval estimate of LHSpline (with ) and GPD (with .95 quantile). Each black (blue) horizontal line is the 90 % CI using the delta (profile likelihood) method for simulated data. The darkgreen horizontal lines are the 90 % CI of the smoothing parameter adjusted LHSpline fit using the conditional simulation approach described in Sec 2.4. The vertical red line denotes the true 50-year return level. The gray horizontal dashed lines are used to divide 10 different Monte Carlo draws.

We summarize in Table 1 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 for 25-, 50-, and 100-year return levels, respectively.

Method Delta Proflik Cond-sim
25-yr   RL ECP 0.63 (46.9) 0.80 (66.9) 0.87 (111.1)
50-yr   RL ECP 0.74 (77.8) 0.84 (93.4) 0.87 (177.9)
100-yr RL ECP 0.77 (107.8) 0.82 (128.9) 0.88 (295.4)
Table 1: The empirical coverage probability (ECP) and (90%) CI width in parentheses for each method for 25-, 50-, and 100-year return levels, respectively.

We also assess the estimation performance of LHSpline in the bulk distribution. Fig 8 shows that the LHSpline performs much better than the gamma fit, a widely used distribution for modeling precipitation (Katz, 1977; Wilks, 2011), and performs nearly as well as the “oracle” (EGPD) approach.

Figure 8: Boxplots of quantile (from 0.001 to 0.999) estimates. The true values are represented by horizontal red lines.

4 Applications

4.1 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 (Risser and Wehner, 2017). In this section, we apply the LHSpline to the daily precipitation measurements for Houston Hobby airport (see Fig. 1 for the histogram and Fig. 9 for the time series) from the Global Historical Climatology Network (GHCN) (Menne et al., 2012).

4.2 Houston Hobby Rainfall Data Analysis

We fit a LHSpline with equally spaced bins to the full range of non-zero precipitation for Hobby Airport ( 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 empirical quantile of nonzero daily precipitation, see diagnostic plot in Fig. 10).

Figure 9: The time series plot of Houston Hobby Airport daily precipitation from 1949 to 2017. Three out of the five largest daily precipitation values were observed during 26th to 28th August, 2017.
Figure 10: Left: The mean residual life plot of the Hobby rainfall data. One should choose the smallest threshold, , such that the mean excess, as a function of , behaves linearly. Right: The parameter stability plot for the shape parameter of GPD of the Hobby rainfall. One should choose the such that the shape parameter estimates stabilize for . The red vertical dotted line in each panel is the chosen threshold in this study.

An important practical issue when applying the LHSpline is that the discretization of rainfall measurement ( of an inch) introduces an artifact in the finely binned histogram in the log scale (see Fig. 11, 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 ( and ) to alleviate this “discreteness” effect.

Figure 11: The density estimate for log scale (left) and original scale (right) of the Hobby airport daily precipitation. The purple dashed (solid) lines are the estimates obtained by LHSpline with (without) adjustment, green dot-dash by LHSpline with lower bound 1 (), golden longdash by LHSpline with lower bound 2 (), red dotted by GPD, and blue dotted by EGPD The red vertical lines indicate the threshold in the log scale (left) and the original scale (right). The “discreteness” of the histogram in the log scale is due to the precipitation measurement precision.

As has been demonstrated in Fig. 3 and Fig. 4 in Sec. 3, a good visual agreement in log-density (density) does not necessarily imply that it actually provides a good return level estimate. Fig. 12 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. 13 (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, , and )) 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. 13 (right panel) also indicates that the LHSpline might provide a better fit than the EGPD to the Hobby Airport data.

Figure 12: The estimated log densities in the log scale (left) and return levels (right). Light red lines are GPD estimates with different thresholds ranging from to empirical quantile. The black step function is the empirical return levels.
Figure 13: The quantile-quantile plot of the Hobby daily precipitation data. Left: A comparison between LHSplines and GPD for the upper 5% () 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:

Method POT LHSpline LHSpline1 LHSpline2 EGPD1 EGPD2
Estimate (years) 30.5 34.8 64.0 98.0 92.5 154.5
90% CI Lower limit 17.0 15.2 19.1 22.0 71.5 117.6
90% CI Upper limit 73.3 73.2 172.4 345.7 457.6 515.6
Table 2: The estimated return period of the 26th August, 2017 daily total precipitation observed at Hobby Airport using POT, LH-Spline, and EGPD.

The much shorter return period estimate obtained from the POT method is due to the overestimation of the upper tail (see Fig. 13) 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. 2).

5 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 (Gilleland and Katz, 2006). 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 or and are largely confined in a bounded region (i.e. , 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 (e.g. Mearns et al., 2009; Wang and Kotamarthi, 2015) 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.


This work was conducted as part of the Research 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.


  • AghaKouchak et al. (2012) A. AghaKouchak, D. Easterling, K. Hsu, S. Schubert, and S. Sorooshian. Extremes in a changing climate detection, analysis and uncertainty. Water Science and Technology Library. Springer, Dordrecht ; New York, 2012. ISBN 1283740990.
  • Alexander et al. (2006) L. Alexander, X. Zhang, T. Peterson, J. Caesar, B. Gleason, A. Klein Tank, M. Haylock, D. Collins, B. Trewin, F. Rahimzadeh, et al. Global observed changes in daily climate extremes of temperature and precipitation. Journal of Geophysical Research: Atmospheres, 111(D5):D05109, 2006.
  • Balkema and De Haan (1974) A. A. Balkema and L. De Haan. Residual life time at great age. The Annals of Probability, 2(5):792–804, October 1974. ISSN 00911798.
  • Boneva et al. (1971) L. I. Boneva, D. Kendall, and I. Stefanov. Spline transformations: Three new diagnostic aids for the statistical data-analyst. Journal of the Royal Statistical Society. Series B (Methodological), 33(1):1–71, January 1971. ISSN 00359246.
  • Brown et al. (2004) L. D. Brown, A. V. Carter, M. G. Low, and C.-H. Zhang.

    Equivalence theory for density estimation, Poisson processes and Gaussian white noise with drift.

    The Annals of Statistics, 32(5):2074–2097, October 2004. ISSN 00905364.
  • Carreau and Bengio (2009) J. Carreau and Y. Bengio. A hybrid pareto model for asymmetric fat-tailed data: the univariate case. Extremes, 12(1):53–76, March 2009. ISSN 1386-1999.
  • Coles (2001) S. Coles. An introduction to statistical modeling of extreme values. Springer series in statistics. Springer, London ; New York, 2001. ISBN 1852334592.
  • Cooley et al. (2007) D. Cooley, D. Nychka, and P. Naveau. Bayesian spatial modeling of extreme precipitation return levels. Journal of the American Statistical Association, 102(479):824–840, September 2007. ISSN 0162-1459.
  • Davison (1984) A. C. Davison. Modelling excesses over high thresholds, with an application. In de Oliveira J.T., editor, Statistical extremes and applications, pages 461–482. Springer, Dordrecht, 1984. ISBN 978-90-481-8401-9.
  • Davison and Smith (1990) A. C. Davison and R. L. Smith. Models for exceedances over high thresholds (with discussion). Journal of the Royal Statistical Society, Series B, Methodological, 52(3), January 1990. ISSN 0035-9246.
  • Dey and Yan (2016) D. K. Dey and J. Yan. Extreme value modeling and risk analysis: methods and applications. Chapman and Hall/CRC, Boca Raton, FL, 2016. ISBN 1498701299.
  • Dupuis (1999) D. Dupuis. Exceedances over high thresholds: A guide to threshold selection. Extremes, 1(3):251–261, January 1999. ISSN 1386-1999.
  • Easterling et al. (2000) D. R. Easterling, G. A. Meehl, C. Parmesan, S. A. Changnon, T. R. Karl, and L. O. Mearns. Climate extremes: observations, modeling, and impacts. Science, 289(5487):2068–2074, 2000.
  • Efron (1979) B. Efron. Bootstrap methods: another look at the jackknife. Annals of Statistics, 7(1):1–26, 1979.
  • Efron and Tibshirani (1996) B. Efron and R. Tibshirani. Using specially designed exponential families for density estimation. The Annals of Statistics, 24(6):2431–2461, December 1996. ISSN 00905364.
  • Eilers and Marx (1996) P. H. Eilers and B. D. Marx. Flexible smoothing with -splines and penalties. Statistical Science, 11(2):89–102, May 1996. ISSN 08834237.
  • Embrechts et al. (1997) P. Embrechts, C. Klüppelberg, and T. Mikosch. Modelling extremal events for insurance and finance. Applications of mathematics ; 33. Springer, New York, 1997. ISBN 3540609318.
  • Fisher and Tippett (1928) R. A. Fisher and L. H. C. Tippett. Limiting forms of the frequency distribution of the largest or smallest member of a sample. Mathematical Proceedings of the Cambridge Philosophical Society, 24(2):180–190, April 1928. ISSN 0305-0041.
  • Frigessi et al. (2002) A. Frigessi, O. Haug, and H. Rue. A dynamic mixture model for unsupervised tail estimation without threshold selection. Extremes, 5(3):219–235, September 2002. ISSN 1386-1999.
  • Gilleland and Katz (2006) E. Gilleland and R. W. Katz. Analyzing seasonal to interannual extreme weather and climate variability with the extremes toolkit. In 18th Conference on Climate Variability and Change, 86th American Meteorological Society (AMS) Annual Meeting, volume 29. Citeseer, 2006.
  • Gilleland and Katz (2016) E. Gilleland and R. W. Katz. extRemes 2.0: An extreme value analysis package in R. Journal of Statistical Software, 72(8):1–39, 2016. doi: 10.18637/jss.v072.i08.
  • Gnedenko (1943) B. Gnedenko. Sur la distribution limite du terme maximum d’une serie aleatoire. Annals of Mathematics, 44(3):423–453, July 1943. ISSN 0003486X.
  • Good and Gaskins (1971) I. J. Good and R. A. Gaskins. Nonparametric roughness penalties for probability densities. Biometrika, 58(2):255–277, August 1971. ISSN 00063444.
  • Green and Silverman (1994) P. J. Green and B. W. Silverman. Nonparametric regression and generalized linear models : a roughness penalty approach. Monographs on statistics and applied probability (Series) ; 58. Chapman and Hall, London ; New York, 1st ed. edition, 1994. ISBN 0412300400.
  • Gu (2013) C. Gu. Smoothing spline ANOVA models. Springer series in statistics, 297. Springer, New York, 2nd ed. edition, 2013. ISBN 1299337511.
  • Gu and Qiu (1993) C. Gu and C. Qiu. Smoothing spline density estimation: Theory. The Annals of Statistics, pages 217–234, 1993.
  • Gumbel (1958) E. J. Gumbel. Statistics of extremes. Columbia University Press, New York, 1958. ISBN 0231021909.
  • Huang et al. (2016) W. K. Huang, M. L. Stein, D. J. McInerney, S. Sun, and E. J. Moyer. Estimating changes in temperature extremes from millennial-scale climate simulations using generalized extreme value (GEV) distributions. Advances in Statistical Climatology, Meteorology and Oceanography, 2(1), July 2016. ISSN 2364-3587.
  • Jenkinson (1955) A. F. Jenkinson. The frequency distribution of the annual maximum (or minimum) values of meteorological elements. Quarterly Journal of the Royal Meteorological Society, 81(348):158–171, April 1955. ISSN 0035-9009.
  • Katz (1977) R. W. Katz. Precipitation as a chain-dependent process. Journal of Applied Meteorology, 16(7):671–676, 1977.
  • Katz et al. (2002) R. W. Katz, M. B. Parlange, and P. Naveau. Statistics of extremes in hydrology. Advances in Water Resources, 25(8):1287–1304, 2002. ISSN 0309-1708.
  • Kharin et al. (2007) V. V. Kharin, F. W. Zwiers, X. Zhang, and G. C. Hegerl. Changes in temperature and precipitation extremes in the IPCC ensemble of global coupled model simulations. Journal of Climate, 20(8):1419–1444, 2007.
  • Kleiber et al. (2012) W. Kleiber, R. W. Katz, and B. Rajagopalan. Daily spatiotemporal precipitation simulation using latent and transformed gaussian processes. Water Resources Research, 48(1), 2012.
  • Kooperberg and Stone (1991) C. Kooperberg and C. J. Stone. A study of logspline density estimation. Computational Statistics and Data Analysis, 12(3):327–347, 1991. ISSN 0167-9473.
  • Leonard (1978) T. Leonard. Density estimation, stochastic processes and prior information. Journal of the Royal Statistical Society. Series B (Methodological), 40(2):113–146, January 1978. ISSN 00359246.
  • Li et al. (2012) C. Li, V. P. Singh, and A. K. Mishra.

    Simulation of the entire range of daily precipitation using a hybrid probability distribution.

    Water resources research, 48(3), 2012.
  • Lindsey (1974a) J. K. Lindsey. Comparison of probability distributions. Journal of the Royal Statistical Society. Series B (Methodological), 36(1):38–47, January 1974a. ISSN 00359246.
  • Lindsey (1974b) J. K. Lindsey. Construction and comparison of statistical models. Journal of the Royal Statistical Society. Series B (Methodological), 36(3):418–425, January 1974b. ISSN 00359246.
  • Martins and Stedinger (2000) E. S. Martins and J. R. Stedinger. Generalized maximum-likelihood generalized extreme-value quantile estimators for hydrologic data. Water Resources Research, 36(3):737–744, 2000.
  • Martins and Stedinger (2001) E. S. Martins and J. R. Stedinger. Generalized maximum likelihood pareto-poisson estimators for partial duration series. Water Resources Research, 37(10):2551–2557, 2001.
  • Mearns et al. (2009) L. O. Mearns, W. Gutowski, R. Jones, R. Leung, S. McGinnis, A. Nunes, and Y. Qian. A regional climate change assessment program for North America. Eos, Transactions American Geophysical Union, 90(36):311–311, September 2009. ISSN 0096-3941.
  • Menne et al. (2012) M. J. Menne, I. Durre, R. S. Vose, B. E. Gleason, and T. G. Houston. An overview of the global historical climatology network-daily database. Journal of Atmospheric and Oceanic Technology, 29(7):897–910, 2012.
  • Naveau et al. (2016) P. Naveau, R. Huser, P. Ribereau, and A. Hannart. Modeling jointly low, moderate, and heavy rainfall intensities without a threshold selection. Water Resources Research, 52(4):2753–2769, April 2016. ISSN 0043-1397.
  • Papastathopoulos and Tawn (2013) I. Papastathopoulos and J. A. Tawn. Extended generalised Pareto models for tail estimation. Journal of Statistical Planning and Inference, 143(1):131–143, January 2013. ISSN 0378-3758.
  • Pickands (1975) J. Pickands. Statistical inference using extreme order statistics. The Annals of Statistics, 3(1):119–131, January 1975. ISSN 00905364.
  • Risser and Wehner (2017) M. D. Risser and M. F. Wehner. Attributable human-induced changes in the likelihood and magnitude of the observed extreme precipitation during hurricane Harvey. Geophysical Research Letters, 44(24):12,457–12,464, December 2017. ISSN 0094-8276.
  • Scarrott and MacDonald (2012) C. Scarrott and A. MacDonald. A review of extreme value threshold estimation and uncertainty quantification. REVSTAT–Statistical Journal, 10(1), March 2012. ISSN 1645-6726.
  • Shaby and Reich (2012) B. A. Shaby and B. J. Reich. Bayesian spatial extreme value analysis to assess the changing risk of concurrent high temperatures across large portions of European cropland. Environmetrics, 23(8):638–648, 2012.
  • Silverman (1982) B. W. Silverman. On the estimation of a probability density function by the maximum penalized likelihood method. The Annals of Statistics, 10(3):795–810, September 1982. ISSN 00905364.
  • Silverman (1986) B. W. Silverman. Density estimation for statistics and data analysis. Monographs on statistics and applied probability (Series). Chapman and Hall, London; New York, 1986. ISBN 0412246201.
  • Smith (1984) R. L. Smith. Threshold methods for sample extremes. In de Oliveira J.T., editor, Statistical Extremes and Applications, pages 621–638. Springer, Dordrecht, 1984. ISBN 978-90-481-8401-9.
  • Smith (1986) R. L. Smith. Extreme value theory based on the r largest annual events. Journal of Hydrology, 86(1-2):27–43, 1986.
  • Tancredi et al. (2006) A. Tancredi, C. Anderson, and A. O’Hagan. Accounting for threshold uncertainty in extreme value estimation. Extremes, 9(2):87–106, June 2006. ISSN 1386-1999.
  • Tapia (1978) R. A. Tapia. Nonparametric probability density estimation. Johns Hopkins series in the mathematical sciences; 1. Johns Hopkins University Press, Baltimore, 1978. ISBN 0801820316.
  • Tebaldi et al. (2006) C. Tebaldi, K. Hayhoe, J. M. Arblaster, and G. A. Meehl. Going to the extremes. Climatic change, 79(3-4):185–211, 2006.
  • Sullivan (1988) F. Sullivan. Fast computation of fully automated log-density and log-hazard estimators. SIAM Journal on Scientific and Statistical Computing, 9(2):363–379, 1988. ISSN 0196-5204.
  • Tsay (2005) R. S. Tsay. Analysis of financial time series, volume 543. John Wiley & Sons, 2005.
  • Wadsworth and Tawn (2012) J. L. Wadsworth and J. A. Tawn. Likelihood-based procedures for threshold diagnostics and uncertainty in extreme value modelling. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 74(3):543–567, June 2012. ISSN 1369-7412.
  • Wahba (1990) G. Wahba. Spline models for observational data. CBMS-NSF Regional Conference series in applied mathematics; 59. Society for Industrial and Applied Mathematics (SIAM, 3600 Market Street, Floor 6, Philadelphia, PA 19104), Philadelphia, PA., 1990. ISBN 0898712440.
  • Wand et al. (1991) M. P. Wand, J. S. Marron, and D. Ruppert. Transformations in density estimation. 86(414):343–353, June 1991. ISSN 0162-1459.
  • Wang and Kotamarthi (2015) J. Wang and V. R. Kotamarthi. High-resolution dynamically downscaled projections of precipitation in the mid and late 21st century over North America. Earth’s Future, 3(7):268–288, 2015. ISSN 2328-4277. 2015EF000304.
  • Weissman (1978) I. Weissman. Estimation of parameters and large quantiles based on the k largest observations. Journal of the American Statistical Association, 73(364):812–815, 1978.
  • Wilks (2011) D. S. Wilks. Statistical methods in the atmospheric sciences. International geophysics series ; v. 100. Academic Press, Oxford ; Waltham, MA, 3rd ed.. edition, 2011. ISBN 0123850223.
  • Zwiers and Kharin (1998) F. W. Zwiers and V. V. Kharin. Changes in the extremes of the climate simulated by CCC GCM2 under doubling. Journal of Climate, 11(9):2200–2222, 1998. ISSN 0894-8755.