Abstract

A hierarchical Bayesian framework is developed to identify multiple abrupt regime shifts in an extreme event series. Specifically, extreme events are modeled as a Poisson process with a gamma-distributed rate. Multiple candidate hypotheses are considered, under each of which there presumably exist a certain number of abrupt shifts of the rate. A Bayesian network involving three layers—data, parameter, and hypothesis—is formulated. A reversible jump Markov chain Monte Carlo (RJMCMC) algorithm is developed to calculate posterior probability for each hypothesis as well its associated within-hypothesis parameters. Based on the proposed RJMCMC algorithm, a simulated example is designed to illustrate the effectiveness of the method. Subsequently, the algorithm is applied to three real, rare event time series: the annual typhoon counts over the western North Pacific (WNP), the annual extreme heavy rainfall event counts at the Honolulu airport, and the annual heat wave frequency in the Chicago area. Results indicate that the typhoon activity over the WNP is very likely to have undergone a decadal variation, with two change points occurring around 1972 and 1989 characterized by the active 1960–71 epoch, the inactive 1972–88 epoch, and the moderately active 1989–2006 epoch. For the extreme rainfall case, only one shift around 1970 is found and heavy rainfall frequency has remained stationary since then. There is no evidence that the rate of the annual heat wave counts in the Chicago area has had any abrupt change during the past 50 years.

1. Introduction

Extreme events are commonly perceived as events that depart pronouncedly from the mean condition. Examples of extreme events include winter storms, hurricanes, heavy rainfall and associated floods, summer heat waves, and freezing spells, to name a few. Because of their potential to cause damage, their occurrences are naturally a matter of concern to the society or ecosystems. As climate is changing, it is of great interest to explore whether there is any pronounced change in the frequency and intensity of extreme events from the past history. Abrupt changes are characteristics of the climate systems. From local, regional, or global perspectives, one of the scientific problems pertinent to our understanding of extreme events is when the abrupt shift occurred and what is the likelihood of its occurrence? In this study, we will address this issue using an advanced statistical modeling approach and apply it to three different kinds of extreme events, namely, tropical cyclones over the western North Pacific; heavy rainfall in Honolulu, Hawaii; and heat waves in Chicago.

Tropical cyclones are one of the most destructive natural catastrophes that recur frequently on the eastern Asian coasts and the Pacific islands. Strong winds, torrential rain, and coastal storm surges often lead to loss of life and enormous property damage. For example, typhoons Rusa in 2003 and Maemi in 2004 inflicted major damage to Korea with a combined damage of more than $10 billion. Heavy rainfall (not necessarily associated with tropical cyclones) is another common extreme event. In Hawaii, heavy rainfall events have often resulted in damage to various sectors with considerable economic loss (Chu et al. 2009). As the earth becomes warmer, it is expected that there are increasing episodes of warm spells. It is well known that increasing mortality results from hyperthermia and cardiovascular disease during heat waves. In Chicago, a total of more than 700 deaths in 1995 were identified as heat related (http://www.noaawatch.gov/themes/heat.php). Given the observed evidence of warmer climate and the vulnerability of our society to extreme events, it is natural to ask whether there is any abrupt shift in the frequency of their occurrences. Have typhoons, heavy rainfall, and summer heat waves become more (or less) active in recent years?

a. Background

The occurrence of rare events is commonly assumed to follow a Poisson process (e.g., Elsner and Bossak 2001; Briggs 2008; Khaliq et al. 2007). For the case of tropical cyclones, Chu and Zhao (2004) developed a three-level hierarchical Bayesian changepoint model, involving data, parameter, and hypothesis layers. In their study, seasonal tropical cyclone counts are represented by a hierarchical Poisson process with gamma-distributed rates. When dealing with this hierarchical model in Bayesian context, it is often very difficult to analytically evaluate complex integrals related to posterior distributions. One efficient way to overcome such difficulty is through the use of the Markov chain Monte Carlo (MCMC) method (e.g., Gamerman and Lopes 2006; Gilks et al. 1996; Gelman et al. 2004). Applications of MCMC to climate research are emerging. For instance, Berliner et al. (2000b) used an MCMC approach to update distribution parameters of their physically based statistical model for predicting the Pacific sea surface temperatures. Elsner et al. (2004) applied a MCMC approach to detect shifts in the Atlantic major hurricane series. More recently, Zhao and Chu (2006) extended the two-hypothesis hierarchical model used in Chu and Zhao (2004) to a multihypothesis hierarchical model and applied the MCMC approach to detecting shifts in the annual hurricane counts in the eastern North Pacific (ENP). They found that the ENP hurricane series is marked by three regimes: the inactive 1972–81 epoch, the active 1982–98 epoch, and the quiescent 1999–2003 epoch. Zhao and Chu (2006) also devised an empirical approach, called the informative prior estimation (IPE) for setting appropriate prior of Poisson rates, followed by framing it into the Bayesian model competing analysis. As a result, they built a more advanced framework for detecting and quantifying multiple (finite) shifts in a series of hurricane counts.

b. Why RJMCMC?

Although the MCMC imbedded with IPE method (Zhao and Chu 2006) has been shown to be viable for calculating the posterior probability for a multiple hypothesis model, it suffers from a shortcoming. That is, because parameter spaces within different hypotheses are typically different from each other, a simulation has to be run independently for each of the candidate hypotheses. If the hypotheses to be investigated have large dimensions, this strategy is not efficient and is computationally prohibitive. In principle, a standard MCMC algorithm is not appropriate for a model selection problem because different candidate models or hypotheses usually do not share the same parameter sets. To overcome this problem, Green (1995) first introduced the reversible jump Markov chain Monte Carlo (RJMCMC) algorithm as a simultaneous integrative approach to deal with the model selection problem.

The RJMCMC algorithm is so termed as to maintain the detailed balance of an irreducible and aperiodic chain that converges to the correct target measure. With the introduction of an extra pseudorandom variable for each investigated model, the dimension mismatching for a different model is well contained. Since the first introduction of the RJMCMC method, many others (e.g., Godsill 2001) have successfully applied and extended it to a broad range of problems such as variable selection, curve fitting, neural network, etc. One of its main real-world applications is on the detection of multiple structure or regime changes in a physical process (e.g., Rotondi 2002). In this study we shall explore a similar problem. Specifically, based on the general RJMCMC framework, we shall design an algorithm to detect the potential multiple abrupt shifts within an extreme event count series and quantify them with proper probabilities.

This paper is organized as follows. Section 2 discusses the datasets we used in this study. Section 3 outlines the model for describing rare event count series. Section 4 briefly introduces the key concepts of RJMCMC method, based on which a complete Bayesian analysis for detecting multiple abrupt shifts is developed. Analysis results of a simulated time series and three real cases are given in section 5. Section 6 provides the summary and conclusions.

2. Data

Tropical cyclone data are from the Joint Typhoon Warning Center in Honolulu, Hawaii. The study period is 1960–2006. In this study, only supertyphoons are used, which are defined as having maximum sustained 1-min wind speed equal to or greater than 130 kt (66.9 m s−1). Historical daily rainfall records for the Honolulu International Airport are obtained from the National Climatic Data Center Web page. The data period extends from 1949 to 2006. To define extreme rainfall events, we choose a specific threshold (99th percentile) of precipitation days as heavy events (Chu et al. 2009). Here we ignore days of no precipitation or days with precipitation less than 0.01 in. (2.5 mm). For each year, we count the number of events when the threshold is crossed. For heat waves, we adopt the World Meteorological Organization definition. That is, daily maximum temperature with more than five consecutive days exceeds the maximum temperature normal by 5°C, where the normal is the period 1961–90. Here Chicago O’Hare International Airport is chosen because of its availability of long records and the infamous 1995 heat waves that led to large number of heat-related deaths. This dataset extends from 1959 to 2006.

3. Mathematical model of a seasonal rare event count series

A Poisson process is a proper probability model for describing independent, rare event counts (e.g., Wilks 2006). Given the rate parameter λ, the probability mass function (PMF) of h counts occurring in T unit seasons is (Epstein 1985)

 
formula

The Poisson mean is the product of λ and T; so is its variance. Throughout this study, all of our cases are annual count series; therefore we always set T = 1 in (1).

In many real applications adopting the model (1), the rate parameter λ is treated as a random variable instead of a constant to construct a hierarchical framework (e.g., Chu and Zhao 2004; Elsner and Jagger 2004). A functional choice of the prior for λ is its conjugate prior, gamma distribution (Epstein 1985) that is formulated by

 
formula

where the gamma function is defined by . That is, given h counts occurring in T years, if the prior density for λ is gamma distributed with parameters h′ and T ′, the resulted posterior density for λ is also gamma distributed with parameters h + h′ and T + T ′. Referring to Eq. (2), the prior expectation of λ is E[λ] = h′/T ′.

Under the statistical model introduced above, the marginal PMF of h counts occurring in T years when the intensity is gamma distributed with prior parameters h′ and T ′ is a negative binomial distribution (Epstein 1985):

 
formula

where h = 0, 1, … , h′ > 0, T ′ > 0, T > 0, and Pnb(·) stands for negative binomial distribution.

4. RJMCMC approach to detecting multiple change points in a rare event count series

a. MCMC method applied to Bayesian analysis

1) Basic inference problem and MCMC method

Let us assume θ to be the set of the model parameters under a given hypothesis, say H, and h to be the data for analysis. The fundamental Bayesian model is governed by the following mathematical expression:

 
formula

where ∝ means proportional. Here, P(h|θ, H) is the conditional distribution of the data h given the model parameters θ within hypothesis H (or the likelihood of the parameters under hypothesis H) and P(θ|H) is the prior distribution of the parameters under hypothesis H. Note that the term in the denominator, the marginal distribution of the data P(h|H), is not dependent on θ; therefore it plays no role if we only focus on evaluating the posterior probabilities of the parameters under hypothesis H. However, in hypothesis or model selection problems, this quantity becomes critical and it is known as the evidence or the marginal likelihood. The Bayes formula provides the posterior probability P(h|θ, H). To make predictive inference, we rely on the posterior predictive distribution P(ĥ|h, H) = ∫ P(ĥ|θ, H)P(θ|h, H)dθ, where ĥ denotes the prediction.

The MCMC method is one of the efficient algorithms for calculating Bayesian inference. Under a hypothesis H (for simplicity, we will drop the notation of H in the following derivation in this subsection), the general Bayesian analysis method described above essentially involves calculating the posterior expectation

 
formula

where a(θ) can be any function of the within-hypothesis parameter set θ. This expectation, however, is usually very difficult to integrate in most practical models. Alternatively, a numerical way to calculate such an expectation is to use an MCMC algorithm by

 
formula

Here, θ[1], θ[2], … , θ[N] are sampled from a Markov chain (MC) after its convergence, where the MC has P(θ|h) as its stationary distribution. This MCMC approach provides an unbiased estimate for E[a|h] (Ripley 1987).

2) RJMCMC approach to Bayesian hypothesis selection

Regular MCMC methods, such as a Gibbs sampler, have been successfully applied to many real applications (e.g., Zhao and Chu 2006). However, they are not appropriate for a Bayesian model selection problem because different models usually do not share the same model parameter set. Generally speaking, a RJMCMC algorithm involves a Metropolis–Hastings-type algorithm that moves a simulation analysis between different models or hypotheses.

Suppose that candidate hypotheses are enumerable and represented by the set H = {H0, H1, … , HK−1}, we denote the joint set of the parameters under hypothesis Hk by θk, k = 0, 1, … , K − 1. We further introduce a random vector uk, k = 0, 1, … , K − 1, such that for any k′ = 0, 1, … , K − 1, k′ ≠ k, the dimension of {θk, uk} and {θk, uk} can be matched. Then we set θk to be a deterministic function of θk and uk. For the reversible move, we propose a vector uk and set θk to be a deterministic function of θk and uk. Thus, there must be a bijection between {θk, uk} and {θk, uk}, which is defined by (θk, uk) = gk,k(θk, uk). A general RJMCMC algorithm (with given observation dataset h) can be sketched as below (after initialization and assuming that the current accepted hypothesis is Hk):

  1. Propose a visit to hypothesis Hk, k′ ≠ k, with probability J(HkHk).

  2. Sample uk from a proposal density Q(uk|Hk, Hk).

  3. Set (θk, uk) = gk,k(θk, uk).

  4. Calculate the odds ratio r and accept Hk as the hypothesis in the next iteration with the probability “min(1, r),” where 
    formula
    If Hk is rejected, we maintain the current hypothesis Hk.
  5. Return to the step 1 until the required number of iterations is reached.

After a burn-in period, the number of times of a hypothesis Hk accepted in the simulation (after step 4) divided by the total number of iterations will be a good estimation for the posterior probability of Hk, P(Hk|h). Also, the samples from each iteration within the hypothesis Hk will be equivalently drawn from the posterior joint probability density function (PDF) P(θk|h, Hk).

In the general RJMCMC algorithm, the choice of an appropriate proposal probability function Q(uk|θk, Hk, Hk) is critical for the efficiency of the algorithm. Many papers have been addressing this issue (e.g., Rotondi 2002). In the following subsections, we shall apply this general RJMCMC framework to the problem on detection of multiple abrupt shifts within a seasonal rare event count series.

b. Hierarchical Bayesian model for extreme event series

Let us assume that there are totally K possible candidate hypotheses {H0, H1, … , HK−1}. Under each hypothesis Hk, k = 0, 1, … , K − 1, there presumably exist exactly k change points in this period. The following derivations are based on the mathematical model described in section 3 and we denote the parameter set under hypothesis Hk by θk. Note that the annual typhoon counts, h = [h1, h2, … , hn]′, are assumed to be a series of independent random variables. The model of hypothesis Hk, k = 0, 1, … , K − 1, is postulated below.

Hypothesis Hk: “A k changes of the rate” of the annual typhoon series:

  •  Assuming that τ = [τk1, τk2, … , τkk]′ denotes the vector of all k change points occurred in the time series h under hypothesis Hk, where the element τk,j represents the first year of the jth epoch. We have the model 
    formula

In (5), τk0 = 1 and τk,k+1 = n + 1 for all k; and the prior knowledge of the hyperparameter set {hkj, T ′kj|j = 1, 2, … , k + 1; and k = 0, 1, … , K − 1} must be given a priori. Note that, under hypothesis Hk, there are k + 1 epochs separated by k change points τkj, j = 1, … , k. The jth change point is defined as the first year of the (j + 1)th epoch.

c. Prior specification

For the prior distribution of the parameter set θk under hypothesis Hk defined in (5), the changepoint parameters are all uniformly (in order) distributed. It is, however, not appropriate to use the noninformative prior for the Poisson rate parameters. In a hypothesis selection problem, a flat noninformative prior usually does not work well because it would almost always favor the simplest hypothesis because of the extremely huge normalization term for each rate parameter. As thoroughly discussed in MacKay (2003), to fit a set of data by using two different models, the posterior probability of the complicated model is penalized by a stronger Occam’s factor, which conceptually is the ratio of its parameters’ posterior and prior widths. Apparently, in model (5), a noninformative prior for the rate will lead to infinite prior width or, in other words, an infinitely small Occam’s factor. To make a sound model or hypothesis selection, a reasonable informative prior for the model parameter is therefore needed. For instance, Zhao and Chu (2006) proposed an approach named IPE to empirically estimate the hyperparameters hkj and T ′kj within each candidate hypothesis. However, IPE involves an extra pilot run for each investigated hypothesis. The accuracy of IPE is also dependent on how well a gamma-Poisson hierarchical model can explain the target time series. Throughout this study, we instead assume that all hkj and T ′kj in model (5) are equal to constant h′ and T ′, respectively. Berliner et al. (2000a) discussed the possibility of using additional information, independent from the time series of interest, by assigning prior distributions. A procedure to fit the hyperprior parameters for a rare event count series is briefly introduced as follows.

With time series h = [h1, h2, … , hn]′, we run L independent iterations. Within the jth iteration, 1 ≤ jL, first we randomly pick two different points from 1 to n, say, k0 and k1(k0 < k1). Then we calculate the sample mean of this batch of samples {hi, k0ik1}, obtaining a realization of the Poisson rate of this iteration, . In the end, we obtain a set of samples, {λ[j], 1 ≤ jL}. Empirically, we assume this Poisson rate is gamma distributed with parameters h′ and T ′. Via a moment estimation approach, an approximation of the hyperparameters h′ and T ′ is given by

 
formula

In this study, after applying the above procedure to several real-world rare event time series (such as seasonal typhoon counts, seasonal heat wave counts, seasonal heavy rainfall count, etc.), we suggest that, for T ′, any number around 10–25 should be a reasonable setting. Throughout this paper, we fix T ′ equal to 18. Analogously, we determine a range for h′, and so we fix h′ = hT ′, where h is the average value of the whole time series h.

d. Extreme event count series analysis by using RJMCMC

As discussed in section 4a, the key for efficiently applying a RJMCMC algorithm is to build an appropriate proposal transition function. For this purpose, we further assume that the hypothesis space defined in model (5) follows a nested structure. In details, referring to step 2 of the general RJMCMC algorithm given in section 4a, to move from the hypothesis Hk to the hypothesis Hk+1, we leave all change points and most rates under Hk untouched and only introduce one new change point, say, τ. Without losing generality, we assume τkj < τ < τk,j+1, 0 ≤ jk. We then propose two new rate parameters for the period [τkj, τ − 1] and [τ, τk,j+1 − 1], respectively, under the hypothesis Hk+1, while discarding the old rate parameter λk,j+1. Conversely, to move from the hypothesis Hk+1 to the hypothesis Hk, we randomly pick one change point under Hk+1, say, τk+1,j, 1 ≤ jk + 1, and then merge the two phases separated by this change point, while introducing a new rate λk,j for the new phase [τk+1,j−1, τk+1,j+1 − 1] under the hypothesis Hk. Simultaneously, we throw away the two old rate parameters and the chosen changepoint parameter. Apparently, with this nested structure model, any two adjacent hypotheses in (5) share most of their parameters.

For the hypothesis transition function J(HkHk), we apply the restriction that any hypothesis can only move to its adjacent hypothesis. A natural set for this is

 
formula

As for the proposal probability Q in (4), generally speaking, it is always preferred that it is close to the posterior probability function of the proposed parameter set under the new hypothesis. A precise estimation of the posterior probability function often involves a pilot run. Nonetheless, with the nested model suggested above, a pilot period is not necessary as an analytical form is available.

First, let us focus on the move from Hk to Hk+1. As we adopt a nested model, hypothesis Hk+1 will retain all k change points under hypothesis Hk. Therefore, this move can be viewed as “birth” of a new change point. Obviously, this “birth” move has the constraint that the new change point cannot take any of the “existent” change points under Hk or the first point of the series, which implies that it has only nk − 1 possible choices. Hence, the noninformative prior for this new change point is P(τ) = 1/(nk − 1). Without losing generosity, we assume that this newly introduced change point τ is within the range [τkj, τk,j+1 − 1], where j can be any integer within the rage [0, k]. In formula, it is τkj < τ < τk,j+1. We then denote the two new rate parameters to be λ1 and λ2 for the two new phases under Hk+1, [τkj, τ − 1], and [τ, τk,j+1 − 1], respectively. Given the new change point τ, the posterior PDFs for λ1 and λ2 under the competing hypothesis Hk+1 are their conjugate:

 
formula

As discussed for the general RJMCMC framework, finding a probability function close to P(τ|θk, h) under Hk+1 is therefore crucial for building an efficient proposal function Qk,k+1. Motivated by the conditional posterior PMF of the new change point τ under Hk+1 (appendix A),

 
formula

we suggest a proposal function P(τ|θk) by substituting λ1 and λ2 in the formula above with their conditional expectation obtained from (8a). It yields

 
formula

We therefore have the proposal probability function from Hk to Hk+1:

 
formula

In (8c), the probability functions P(λ1|θk, τ), P(λ2|θk, τ), and P(τ|θk) are calculated via (8a) and (8b).

Referring to (4), we also need to calculate the proposal probability for the reverse move from hypothesis Hk+1 to hypothesis Hk, in which we need to propose a new rate, denoted by λ, for the period [τkj, τk,j+1 − 1] under Hk or, equivalently, [τk+1,j, τk+1,j+2 − 1] under Hk+1. The posterior PDF for λ under Hk is its conjugate:

 
formula

We then can draw λ from (9).

Referring to step 3 in the general RJMCMC algorithm, the bijection function is the identity mapping, which implies that Jacobian determinant term is equal to 1. The likelihood term can be found in appendix B. Thus, after the overlapped parts between two hypotheses are canceled out with each other, the ratio in (4) is simplified to (assuming τkj < τ < τk,j+1)

 
formula

In (10a), P(λ), P(λ1), and P(λ2) are all gamma distributed with hyperparameters h′ and T ′, where h′ and T ′ are obtained through the procedure described in section 4c [Eq. (6)]; P(τ) = 1/(nk − 1); hypothesis transition probability J is defined in (7); proposal functions Qk+1,k(λ|θk+1) and Qk,k+1(λ1, λ2, τ|θk) are given in (9) and (8c), respectively. The derivation of the likelihood within each hypothesis is found in appendix B.

For the reverse move from hypothesis Hk+1 to hypothesis Hk, we simply execute the opposite transition from Hk to Hk+1. That is, we shall merge two adjacent periods under Hk+1 into one period with this move, which can be viewed as “death” of a change point under hypothesis Hk+1. Without losing generosity, assume the “dead” change point to be τk+1,j for any applicable j. We thus first sample a new rate λ = λk,j via (9) for the new merged phase, [τk+1,j−1, τk+1,j+1 − 1] (or equivalently [τkj, τk,j+1 − 1] under hypothesis Hk). For the switch back, we sample λ1, λ2, and τ from (8a) and (8b), respectively, within this newly merged period. As a “death” move is exactly the reverse of its “birth” counterpart, the acceptance ratio for this move simply takes the inverse of the form in (10a). That is,

 
formula

Replacing the ratio calculation equation in (4) with (10a) and (10b), the general RJMCMC algorithm provided in section 4a yields the proposed algorithm ad hoc designed for multiple changepoint detection within an extreme event count series. In real practice, for both (10a) and (10b), we further assume a noninformative prior for all candidate hypotheses. That, is all hypotheses have equal probability before analyzing the target series. That is, P(Hk) = 1/K is a constant for k = 0, 1, … , K − 1. The prior term for hypothesis in (10a) and (10b) is thereby canceled out.

5. Results

In this study, we assume that there are at most nine change points, that is, K = 10. For each of the following examples, the length of the burn-in period is 2000 and the number of iterations of the RJMCMC algorithm after the burn-in period is chosen as 10 000. To test the convergence of RJMCMC simulation chains, we apply the diagnostic software Bayesian Output Analysis (BOA; Smith 2007).

a. A simulated example

To test the validity and effectiveness of the proposed RJMCMC method, a simulated example was elaborately designed. We first generated a time series with gamma distribution. There are 500 points of the series in total. The first 150 points were generated from a gamma(3, 1); the next 150 points were generated from a gamma(5, 1); the following 100 points were generated from a gamma(9, 3); the last 100 points were generated from a gamma(3, 0.5). With this series as the seed of the rate, we generate a Poisson-distributed time series. Obviously, with this generative model, there should be four independent epochs, which are separated by three change points with the sample index 151, 301, and 401, respectively. Figure 1a illustrates this simulated series.

Fig. 1.

Plots for the simulated example. (a) Time series. (b) Posterior probability for all the candidate hypotheses. (c) Posterior PMF for the first change point. (d) Posterior PDF of the Poisson rate for the first epoch. (e) Posterior PDF of the Poisson rate for the second epoch. (f) Posterior PMF for the second change point. (g) Posterior PDF of the Poisson rate for the third epoch. (h) Posterior PDF of the Poisson rate for the forth epoch. (i) Posterior PMF for the third change point.

Fig. 1.

Plots for the simulated example. (a) Time series. (b) Posterior probability for all the candidate hypotheses. (c) Posterior PMF for the first change point. (d) Posterior PDF of the Poisson rate for the first epoch. (e) Posterior PDF of the Poisson rate for the second epoch. (f) Posterior PMF for the second change point. (g) Posterior PDF of the Poisson rate for the third epoch. (h) Posterior PDF of the Poisson rate for the forth epoch. (i) Posterior PMF for the third change point.

After applying the RJMCMC method described in the preceding section to this series, the outputs for the posterior probability of 10 candidate hypotheses are demonstrated in Fig. 1b. The probability for hypothesis H3 is as high as 0.82 and the Bayes factor between H3 and all other hypotheses is equal to 4.6, suggesting H3 is the dominant winner. Under hypothesis H3, the posterior PMFs for three change points are illustrated in Figs. 1c, 1f, and 1i, respectively, from which we see the algorithm precisely catches the prescribed abrupt shifts. We also plotted the posterior PDFs of the Poisson rates for each independent epoch in Figs. 1d, 1e, 1g, and 1h, respectively. For this simulated example, very satisfying results were obtained as expected.

b. Three real-world examples

Now let us turn our attention to modeling the actual climate records. This includes three cases: the annual supertyphoon counts over the western North Pacific (WNP); the annual extreme heavy rainfall counts at Honolulu International Airport; and the annual heat wave counts in Chicago’s O’Hare International Airport. To estimate the p value associated with a rate shift under a given hypothesis, we calculate the total number of samples that are against the presumed rate shift trend and divide this number to the total number of iterations under this hypothesis that are accepted by the RJMCMC algorithm (after burn-in period).

1) Annual supertyphoon counts in the WNP

Figure 2a shows the time series of annual typhoon count series over the WNP from 1960 to 2006. After applying the same procedure as described in the first example, the outputs for the posterior probability of 10 candidate hypotheses are plotted in Fig. 2e. The probability for hypothesis H2 is as high as 0.39, which is higher than that for hypothesis H3 (0.24). We therefore choose H2 as the winning hypothesis. Under hypothesis H2, the marginal posterior PMF for the two change points are depicted in Figs. 2c and 2d, respectively, through which one can see that by far the most likely choice for the first change point is 1972 and for the second change point is 1989. Using a method of running medians, Yumoto and Matsuura (2001) analyzed the annual tropical cyclone (TC) records in the western North Pacific. They identified 1961–72 as an active epoch, followed by an inactive epoch of 1973–85 and another active epoch of 1986–94. They attributed this regime shift to the corresponding interdecadal variability in the atmospheric and oceanic conditions in the subtropical western North Pacific. Interestingly, our study also independently confirms a phase shift occurring around 1972. However, besides providing a more rigorous and advanced method in detecting change points, our study reveals the likelihood of the regime shift in any year and the posterior probability of each hypothesis, which are unavailable in Yumoto and Matsuura (2001).

Fig. 2.

Plots for the case of supertyphoon counts series. (a) Time series. (b) Posterior PDFs for the rates of each epoch. (c) Posterior PMF for the first change point. (d) Posterior PMF for the second change point. (e) Posterior probability for all the candidate hypotheses.

Fig. 2.

Plots for the case of supertyphoon counts series. (a) Time series. (b) Posterior PDFs for the rates of each epoch. (c) Posterior PMF for the first change point. (d) Posterior PMF for the second change point. (e) Posterior probability for all the candidate hypotheses.

The posterior PDF for the rates of each epoch, along with the universal prior, are collaboratively demonstrated in Fig. 2b. The average rate prior to 1972 is about 5.67 supertyphoons yr−1 and decreases to about 2.35 supertyphoons yr−1 from 1972 to 1988, then it increases to 5.00 supertyphoons yr−1 thereafter. The p values for the both shifts are all very small (0.054 for the first and 0.059 for the second), which strongly implies the existence of these two change points in this storm time series.

The Bayes factor between H2 and H3 is 1.61, which only shows slight evidence in favor of H2 over H3. However, we observe that, under H3, based on the obtained changepoint samples, besides the two shifts (1972 and 1989) identified under H2, the third most probable change point is 1987, which is very close to 1989. In summary, we suspect that there are two regime shifts that occurred within this series: one was around the early 1970s and the other was around late 1980s.

2) Annual extreme heavy rainfall counts at Honolulu International Airport

Figure 3a shows the time series of the annual extreme heavy rainfall count series observed in Honolulu International Airport from 1950 to 2006. We applied the same procedure described in the last example to this series, and the outputs for the posterior probability of 10 candidate hypotheses are plotted in Fig. 3b. The probability for hypothesis H1 is 0.36, which is higher than that for hypothesis H2 (0.26) or hypothesis H0 (0.19). We therefore chose H1 as the winning hypothesis. Under hypothesis H1, the marginal posterior PMF for the only change point is plotted in Fig. 3d, through which one can see that the decreasing shift more likely happened around 1970.

Fig. 3.

Plots for the case of extreme rainfall counts series. (a) Time series. (b) Posterior probability for all the candidate hypotheses. (c) Posterior PDFs for the rates of both epochs. (d) Posterior PMF for the change point.

Fig. 3.

Plots for the case of extreme rainfall counts series. (a) Time series. (b) Posterior probability for all the candidate hypotheses. (c) Posterior PDFs for the rates of both epochs. (d) Posterior PMF for the change point.

The posterior PDF for the rates of each epoch, along with the universal prior, are all shown in Fig. 3c. The average rate prior to 1970 is about 1.40 times yr−1 and is almost halved to 0.70 times yr−1 since then. The p value for the shift is 0.096, which supports the existence of this change point in this time series. Since the Bayes factor between H1 and H2 is only 1.40, based on Jeffreys’ criterion, this change point needs further investigation to confirm.

3) Annual heat wave counts in Chicago area

Figure 4a shows the time series of the annual extreme heavy rainfall count series observed in O’Hare airport from 1959 to 2006. The output for the posterior probability of 10 candidate hypotheses is plotted in Fig. 4c. The probability for hypothesis H0 is 0.34, which is slightly higher than that for hypothesis H1 (0.31). Under hypothesis H0, the posterior PDF for the rate through the period, along with its prior, are drawn in Fig. 4c. The average rate of this series is about 4.33 times yr−1.

Fig. 4.

Plots for the case of heat wave counts series. (a) Time series. (b) Posterior probability for all the candidate hypotheses. (c) Posterior PDF for this period.

Fig. 4.

Plots for the case of heat wave counts series. (a) Time series. (b) Posterior probability for all the candidate hypotheses. (c) Posterior PDF for this period.

Since the Bayes factor of H0, H1 over all other hypotheses is 1.84, we exclude those higher dimensional hypotheses. With our further analysis on the samples drawn within iterations under H1, we find that the posterior PMF for the only change point under H1 is almost like a uniform distribution. We thus claim that no change point exists within this time series.

6. Summary

As the earth’s climate is changing, the frequency in extreme events is expected to change accordingly. Hence, developing a state-of-the-art method to objectively identify the turnaround of such changes is an initial vital step for a more comprehensive scientific analysis. With the benefit of knowing when a regime shift has occurred, this knowledge would enable one to compare active and inactive epochs of climate states for future diagnostic and modeling studies (e.g., Chu 2002; Deser et al. 2004). Traditionally, rates for seasonal extreme event counts have been modeled in a data-parameter two-layer hierarchical Bayesian framework. In this view, the rates are assumed to be invariant throughout the time. A few studies respectively provided different Bayesian approaches to detecting and quantifying potential abrupt shifts in an extreme event series (e.g., Elsner et al. 2004; Zhao and Chu 2006) by treating the Poisson rate a random variable. However, each of them has some limitations as we have discussed in the introduction section.

A general 3-layer—including data, parameter, and hypothesis—hierarchical Bayesian model with a nested hypothesis space is built in this study. Seasonal extreme event count series is modeled as a Poisson process with a gamma-distributed rate. We consider multiple candidate hypotheses, within each of which there presumably exist a certain number of abrupt shifts of the Poisson rate. Because of the inability of a regular MCMC algorithm for a model selection problem, we resort to its extended version, the RJMCMC algorithm. Following the guideline of the general RJMCMC method, we designed an algorithm to automatically calculate the Bayesian inference of the 3-layer Bayesian model for efficiently solving the hypothesis competition problem.

A simulated example was designed to illustrate the effectiveness and consistency of the proposed RJMCMC method, and satisfying results were obtained. Subsequently, the algorithm was applied to three examples of extreme event series: the annual supertyphoon counts over the WNP, the extremely heavy rainfall counts in Honolulu airport, and the annual heat wave counts in Chicago area. The results indicate that typhoon activity over the WNP is very likely to have undergone a decadal variation with two change points occurring around 1972 and 1989; the average supertyphoon rate is 5.67 yr−1 during the active 1960–71 epoch, drops to 2.35 yr−1 during the inactive 1972–88 epoch, and then goes up to 5.00 yr−1 from 1989 to 2006. The extreme rainfall occurrence frequency in Honolulu airport may have had a decreasing shift in around 1969 and it has remained stationary since then; there is no evidence that the rate of the annual heat wave counts in Chicago area has had any significant variation within the past half century.

One cautionary note is that the mathematical models and computer codes developed here are based on the extreme event series, which are assumed to follow a Poisson process with a nonconstant rate. Because of this assumption, the algorithm presented in this study is not applicable to the nonrare event series. To apply the proposed algorithm in real practice, it is necessary to check the independency of the tested time series. A simple way is to calculate the autocorrelation function and then compare the output correlation coefficients to a threshold value. Furthermore, after running the algorithm, for those slightly significant change points (say Bayes factor smaller than 2), it may be necessary to check other relevant meteorological phenomena to confirm their existence. To make this algorithm be more general, on a side note, besides checking the target time series’ independency, it would be also helpful to test if it has any significant period by using a frequency analysis tool, such as harmonic analysis. How to properly model the series embedded with an oscillation feature, however, is beyond the scope of this study.

Acknowledgments

This study is partially funded by the Pacific Disaster Center. Thanks are also due to the anonymous reviewers for their constructive and insightful suggestions that led to substantial improvements in the manuscript.

REFERENCES

REFERENCES
Berliner
,
L. M.
,
R. A.
Levine
, and
D. J.
Shea
,
2000a
:
Bayesian climate change assessment.
J. Climate
,
13
,
3805
3820
.
Berliner
,
L. M.
,
C. K.
Wilke
, and
N.
Cressie
,
2000b
:
Long-lead prediction of Pacific SSTs via Bayesian dynamic modeling.
J. Climate
,
13
,
3953
3968
.
Briggs
,
W. M.
,
2008
:
On the changes in number and intensity of North Atlantic tropical cyclones.
J. Climate
,
21
,
1387
1402
.
Chu
,
P-S.
,
2002
:
Large-scale circulation features associated with decadal variations of tropical cyclone activity over the central North Pacific.
J. Climate
,
15
,
2678
2689
.
Chu
,
P-S.
, and
X.
Zhao
,
2004
:
Bayesian changepoint analysis of tropical cyclone activity: The central North Pacific case.
J. Climate
,
17
,
4893
4901
.
Chu
,
P-S.
,
X.
Zhao
,
M.
Grubbs
, and
Y.
Ruan
,
2009
:
Extreme rainfall events in the Hawaiian Islands.
J. Appl. Meteor. Climatol.
,
48
,
502
516
.
Deser
,
C. A.
,
S.
Philips
, and
J. W.
Hurrell
,
2004
:
Pacific interdecadal climate variability: Linkages between the tropics and the North Pacific during boreal winter since 1900.
J. Climate
,
17
,
3109
3124
.
Elsner
,
J. B.
, and
B. H.
Bossak
,
2001
:
Bayesian analysis of U.S. hurricane climate.
J. Climate
,
14
,
4341
4350
.
Elsner
,
J. B.
, and
T. H.
Jagger
,
2004
:
A hierarchical Bayesian approach to seasonal hurricane modeling.
J. Climate
,
17
,
2813
2827
.
Elsner
,
J. B.
,
X.
Niu
, and
T. H.
Jagger
,
2004
:
Detecting shifts in hurricane rates using a Markov chain Monte Carlo approach.
J. Climate
,
17
,
2652
2666
.
Epstein
,
E. S.
,
1985
:
Statistical Inference and Prediction in Climatology: A Bayesian Approach.
Meteor. Monogr., No. 42, Amer. Meteor. Soc., 199 pp
.
Gamerman
,
D.
, and
H. F.
Lopes
,
2006
:
Markov Chain Monte Carlo—Stochastic Simulation for Bayesian Inference.
2nd ed. Chapman & Hall/CRC, 323 pp
.
Gelman
,
A.
,
J. B.
Carlin
,
H. S.
Stern
, and
D. B.
Rubin
,
2004
:
Bayesian Data Analysis.
2nd ed. Chapman & Hall/CRC, 668 pp
.
Gilks
,
W. R.
,
S.
Richardson
, and
D. J.
Spiegelhalter
,
1996
:
Markov Chain Monte Carlo in Practice.
Chapman & Hall/CRC, 486 pp
.
Godsill
,
S. J.
,
2001
:
On the relationship between Markov chain Monte Carlo methods for model uncertainty.
J. Comput. Graph. Stat.
,
10
,
1
19
.
Green
,
P.
,
1995
:
Reversible jump Markov chain Monte Carlo computation and Bayesian model determination.
Biometrika
,
82
,
711
732
.
Khaliq
,
M. N.
,
T. B. M. J.
Ouarda
,
A.
St.-Hilaire
, and
P.
Gachon
,
2007
:
Bayesian changepoint analysis of heat spell occurrences in Montreal, Canada.
Int. J. Climatol.
,
27
,
805
818
.
MacKay
,
D. J. C.
,
2003
:
Information Theory, Inference, and Learning Algorithms.
Cambridge University Press, 628 pp
.
Ripley
,
B. D.
,
1987
:
Stochastic Simulation.
John Wiley, 237 pp
.
Rotondi
,
R.
,
2002
:
On the influence of the proposal distributions on a reversible jump MCMC algorithm applied to the detection of multiple changepoints.
Comput. Stat. Data Anal.
,
40
,
633
653
.
Smith
,
B. J.
,
cited
.
2007
:
Bayesian Output Analysis Program (BOA), version 1.1.5 for R.
Department of Biostatistics, School of Public Health, University of Iowa. [Available online at http://www.public-health.uiowa.edu/boa/]
.
Wilks
,
D. S.
,
2006
:
Statistical Methods in the Atmospheric Sciences.
2nd ed. Academic Press, 627 pp
.
Yumoto
,
M.
, and
T.
Matsuura
,
2001
:
Interdecadal variability of tropical cyclone activity in the western North Pacific.
J. Meteor. Soc. Japan
,
79
,
23
35
.
Zhao
,
X.
, and
P-S.
Chu
,
2006
:
Bayesian multiple changepoint analysis of hurricane activity in the eastern North Pacific: A Markov chain Monte Carlo approach.
J. Climate
,
19
,
564
578
.

APPENDIX A

Derivation of P(τ|θk, λ1, λ2, h) within Hypothesis Hk+1

Let h = [h1, h2, … , hn] and presumably τkj < τ < τk,j+1, j ∈ {0, … , k}. With the assumed nested structure, it implies τk,j = τk+1,j < τ = τk+1,j+1 < τk,j+1 = τk+1,j+2. Therefore, the likelihood for the data before the change point τk+1,j and after the change point τk+1,j+2 does not contain any information of τ. With model (5), we have the likelihood function for the dataset:

 
formula

Ignoring the terms that do not contain τ, we have

 
formula

Since P(τ) is assumed to be uniformly distributed, we have

 
formula

In (A2), we simply replace the parameter index within Hk+1 in (A1) by the relative parameter index within Hk.

APPENDIX B

Derivation of Likelihood P(h|θk, Hk) k = 0, 1, … , K − 1

We have

 
formula

Therefore,

 
formula

APPENDIX C

Software Availability

The MATLAB code developed based on this paper, as well the relative help documents, can be downloaded from the corresponding author’s homepage (http://lumahai.soest.hawaii.edu/Hsco/chu-index.html).

Footnotes

Corresponding author address: Dr. Pao-Shin Chu, Department of Meteorology, School of Ocean and Earth Science and Technology, University of Hawaii at Manoa, 2525 Correa Road, Honolulu, HI 96822. Email: chu@hawaii.edu