A Space-Time Model for Seasonal Hurricane Prediction

THOMAS H. JAGGER AND XUFENG NIU
DEPARTMENT OF STATISTICS, THE FLORIDA STATE UNIVERSITY, TALLAHASSEE, FLORIDA
CORRESPONDING AUTHOR ADDRESS: INSIGHTFUL, SEATTLE, WA
TEL: 206-283-8802X367
E-MAIL: TJAGGER@INSIGHTFUL.COM

JAMES B. ELSNER
DEPARTMENT OF GEOGRAPHY, THE FLORIDA STATE UNIVERSITY, TALLAHASSEE, FLORIDA

May 23, 2001

Abstract

1. Introduction

Hurricanes are destructive natural phenomena. On average a half of dozen or so form each year over the warm tropical and subtropical waters of the North Atlantic and track westward for thousands of kilometers. Historically, hurricanes account for a majority of the costliest weather disasters in the United States. They rival earthquakes in destructive potential and loss of life. Despite technological advances in monitoring and prediction, hurricanes retain their potential to cause severe damage and numerous deaths (Arguez and Elsner 2001). During 1998, hurricane Mitch became a grim reminder that hurricanes can quickly kill thousands of people. In the United States, population and demographic shifts toward the coast are making the problem worse as development flourishes in areas prone to hurricane strikes: the warm subtropical shorelines and islands of the Atlantic Ocean and Gulf of Mexico.

Knowledge of their past occurrence, even if it is incomplete, provides clues about their future frequency and intensity that goes beyond what present numerical climate models are capable of. This understanding is important for land-use planning, emergency management, hazard mitigation, (re)insurance applications and, potentially, the long-term weather derivative market. Climatologists have been issuing seasonal hurricane activity forecasts for the North Atlantic using since 1984 (Gray, 1984; Elsner et al., 1996). These forecasts provide total basin estimates of annual counts using linear and generalized linear regression models. Forecasts are issued several months in advance of the season, which runs from June through November, and updated as the season approaches.

Some understanding has been achieved in solving the climate puzzle with regard to the question of where hurricanes are likely to go based on conditions a month or two in advance of the season. Lehmiller et al. (1997) outline the problem and show the potential for specific forecast models based on regional and large-scale climate factors. Elsner et al. (2000) demonstrate a link between the probability of a major hurricane strike on the East Coast and the strength of the North Atlantic oscillation (NAO). Physically it is reasoned that the strength of the NAO is an indication of the position and intensity of the subtropical high, which steers the hurricanes toward the coast. Taking advantage of a unique wind speed data set from hurricane landfalls in the United States, Jagger et al. (2001) develop a model for hurricane probabilities conditioned on climate anomalies, including El Niño and the NAO.

But current seasonal forecast models for the entire North Atlantic basin do not incorporate spatial and temporal information. Thus they fail to provide specific seasonal activity forecasts for different geographic regions of the hurricane basin, which includes the Gulf of Mexico and Caribbean Sea. Statistical models need to be created that combine spatial and temporal correlation in the data to generate regional forecasts. This paper introduces a class of space-time statistical models for count data that can be used for seasonal hurricane prediction. Ultimately the goal is a model that can predict the likely tracks of hurricanes for an entire season.

The paper is organized as follows. Section two describes the model grid and data. Model formulae and justification are presented in section three along with a description of the model predictors. Section four describes increasingly sophisticated models as a way to understand the full model. Final model selection procedures are given in section five. Some issues related to the estimation procedure and a single hindcast are given in section six. The paper ends with a summary in section seven.

2. Model grid and data

To develop the model we divide the North Atlantic basin into a 6^° by 6^° latitude and longitude grid. Grid choice is a compromise between sample size and resolution. As noted later, model coefficients are estimated using a Monte Carlo procedure, and the procedure fails to converge for smaller grid boxes because between grid box correlations are too large.

Data are annual hurricane counts in the grid boxes. A tropical cyclone that records a position at hurricane intensity within the box is counted once. Hurricanes are tropical cyclones with maximum sustained winds reaching 65 kt or greater. A hurricane that loops around and reenters the box is counted as a single hurricane. Hurricane positions and intensities are obtained from the best-track records (Neumann et al., 1999), which are a compilation of the six-hourly information of all tropical cyclones back to 1886. Records are most reliable after 1943, particularly for the weaker, shorter lived cyclones. The period of record for the present study is 1900-93 (94 yr).

We remove grids having mostly land or historically low hurricane activity, leaving grid S with 40 cells as shown in Fig. . Total hurricane occurrences over the 94-yr period are shown for each region. Historically, hurricane activity is most pronounced over the Bahamas extending north and eastward toward Bermuda (see Elsner and Kara 1999). Count values in the grid boxes outside region S are used as boundary conditions for the model. This improves the estimates of model constants over the case where the boundary values must be assumed or estimated. For generating sample forecasts we set the boundary values equal to zero to approximate annual climatology.

Count values in the grid boxes obtained from the best-track data and used in model development are displayed as a series of space-time images in Fig. . Each image represents an 11 by 8 grid (6^° by 6^° latitude/longitude), covering the western half of the North Atlantic for a single season. Gray levels in the images indicate the number of hurricanes whose centers pass into the region during the year. Activity in the grids range from zero (light gray) to four (black). Low latitude hurricane activity dominated the region during 1916 and 1933 whereas high latitude activity was more pronounced during 1963 and 1969. Note that there are no hurricanes in the best-track data for the years 1907 and 1914. These annual maps of seasonal hurricane counts are used as the model response as detailed next.

3. A right-truncated Poisson space-time model

Hurricane frequencies over time and space form a space-time counts process. The dependence structure of this type of data can be modeled by a conditional probability approach (Bartlett, 1968; Whittle, 1963; Besag, 1974; Gilks et al., 1996). Besag (1974) introduce conditionally specified auto-Poisson models for spatial counts data, which link observation of a Poisson process at a given location with those at its spatial neighborhoods. But the auto-Poisson model proposed by Besag (1974) has restrictions on the parameter space making it applicable only to spatial data in which the interaction coefficients are non-positive. Here we consider a class of space-time regression models for hurricane activity based on the right-truncated Poisson distribution

Pr(H=k) = lk k! M å j=0  (lj/j!) ,     k=0, 1, ¼, M, (1)

where M is an upper bound for the annual number H of hurricanes occurring in a given grid box during any given year. The truncated Poisson space-time autoregressive (TRSTAR) model we propose allows the coupling parameters for neighboring sites to take on positive values.

We restrict the specification to a spatially invariant, nearest neighbor, first order autoregressive model, with predictors. The space-time neighborhood consists of five sites: north, south, east, west, and the previous season's activity at the given site as shown in Fig. . Note that using the TPSTAR model on this neighborhood requires three coupling parameters g_0h,g_0v,g_1c for the east-west, the north-south and the temporal lag-one coupling, respectively.

Each grid box (t,ij) represents a 6^° by 6^° region for one year. Let the response values H_t,ij, be the number of distinct hurricanes passing into any portion of the grid (i,j) Î S during year t. The longitude of the location centers are given by -107^° + 6^°·i, for i=1¼11 and the latitude centers by 6^°+6^°·j, for j=1¼8. We can describe the distribution model conditionally as:

"ij Î S  :  Ht,ij | h¶{ (t,ij)} ~ tpois(lt,ij,M) (2)

where h_¶(t,ij) is the vector of five response values, h_t,i-1,j,h_t,i+1,j,h_t,i,j-1,h_t,i,j+1 and h_t-1,i,j, in the neighborhood of the location, (t,i,j), and tpois is the right-truncated Poisson distribution with rate l_t,ij truncated at M. For the neighborhood system shown in Fig. , the TRSTAR model has the form:

log(lt,ij) = aij(t)+ g0h(Ht,i+1,j+Ht,i-1,j)+g0v (Ht,i,j+1+Ht,i,j-1) (3) + g1c (Ht-1,i,j) where aij(t) = aij + p å k=1  ck zk(t)

with c_k the coefficient for the k^th yearly covariate and a_ij the grid ij offset.

We chose M=10 for the grid boxes, since it is larger than the maximum number of hurricanes observed in any grid box per year. The value of M is small enough to observe positive coupling, yet large enough to compare the results of this model with models based on the Poisson distribution.

The model makes use of three predictors types:

1.

neighborhood response values or coupling grid box frequencies,

2.

local offsets, represented by a_ij, and

3.

yearly global covariates denoted by z_k(t).

Predictors of type two represent the regional effects from predictors that are not measured at each location, such as sea surface temperature, location to land mass, and the variation in surface area from region to region. The parameters are given by a_ij. Without additional information the site offset gives an indication of average yearly activity in the associated grid box. As a_ij represents one parameter for every grid box, this can be a considerable number of parameters. Thus, we constrain a_ij to be the sum of two factors, longitude a_i and latitude b_j, where a_ij=a_i+b_j. This reduces the number of parameters associated with offsets to 14 in the present case.

Predictors of type three represent the climate covariates as identified by previous research studies (Jagger et al., 2001). Here we consider the influence of the covariate to be the same across grid boxes, and allow the value of the covariate to vary in intensity from year to year. The number of covariates is restricted by data availability. The climate covariates used in the TPSTAR model include:

Warm & Cold

``Warm" and ``Cold" refer to the state of the El Niño-Southern Oscillation (ENSO) which is the aperiodic warming and cooling of the ocean mixed-layer waters over the equatorial central and eastern tropical Pacific and their associated atmospheric interaction. This factor takes three values in the model, ``Warm", +1 representing an El Niño event; ``Cold", -1, representing a La Niña event or 0 for neither event. Definitions and values are given in Elsner et al. (1998). Shapiro (1987) demonstrates wind speeds at levels above the tropical North Atlantic are correlated to the phase of the ENSO. Since stronger winds in the upper atmosphere are related to increased wind shear over areas of tropical storm development, the presence or absence of El Niño is an indicator of the level of hurricane activity over the North Atlantic.

Dakar

``Dakar" refers to annual rainfall over Dakar, Senegal. Landsea and Gray(1992) show a correlation between North Atlantic hurricane activity and rainfall during the west African rainy season. It is speculated that above normal rainfall in this region is a signal that sea surface temperatures over the tropical North Atlantic are above normal and that therefore the adjacent ocean waters will likely spawn more tropical storms and hurricanes.

Azores

``Azores" refers to annual sea level pressure (SLP) in the Azores given in millibars. Sea level pressure over the Azores is an indication of the strength and position of the subtropical high pressure zone.

Iceland

``Iceland" refers to annual SLP in Iceland given in millibars. Elsner and Kara (1999) note that lower pressures are favorable for hurricane formation, either by reduced wind shear due to relaxed trades or baroclinic enhancement from weak low pressure eddies in the sub-tropics generated by major middle latitude troughs. The difference in SLP between Iceland and the Azores (Iceland - Azores) represents the North Atlantic oscillation. Elsner et al. (2000) show that major hurricanes along the Northeast coast of the United States are more likely under a strong, positive NAO.

4. Comparing models

To better examine hurricane activity as a count process over the basin, we present results from three different models, a time-series only model, a space-time model without instantaneous coupling parameters, and a full model with coupling parameters. The full model given in (3) represents a version of TPSTAR. The first two models are fit to a Poisson generalized linear model (GLM) with dispersion. To fit the full model we try both a maximum Poisson likelihood estimator (MPLE) and a Monte Carlo maximum likelihood estimator (MCMLE).

The MPLE extends the GLM model by adding neighborhood observations as covariates to the models. If the covariate matrix is of full rank, the MPLE using the canonical link function always produces a parameter estimate. If the model is shift invariant with finite range, the parameter estimates of the MPLE are asymptotically consistent with increasing domains (Winkler 1995). However, even in the shift invariant case, the MPLE is not necessarily efficient, and does not provide standard error estimates. We use the MCMLE method for consistent and efficient parameter estimates with consistent parameter covariance estimates as described in Geyer (1994) and used by Wu (1994). The MCMLE method was extended by Jagger (2000) to handle the autoregressive coupling in the TPSTAR model.

The MCMLE method has several problems. For one, it is computationally intensive as it uses a Markov Chain Monte Carlo method for estimating the log likelihood function. For another, the method fails to converge, unless the initial parameter estimates are close to the actual parameter values. Though the MPLE is biased in our case, the convergence problem is somewhat alleviated by using the MPLE for initial parameter estimates.

For the time-series model, we analyze annual hurricane counts over the entire region S. The time series for 95 years, H_t, t Î 1889 ¼1993 is shown in Fig. . H_t represents the number of separate hurricanes passing through any portion of the region in year t. Although not all predictors are significant in the time-series formulation, we compare the parameter values and their errors to values obtained after dividing the region into 6^° by 6^° grid boxes.

Table  shows the parameter estimates from the time series model. The formulation is a Poisson GLM with dispersion. A Poisson GLM is used on U.S. hurricane activity in Elsner and Bossak (2001). We note that the model is somewhat under dispersed. The dispersion is estimated from the Pearson statistic å_t=1ⁿ[(H_t-f_t)/SE(f_t)]² where f_t is the fitted value for H_t and SE(f_t) is the standard error of f_t. Dakar rainfall, Azores pressure, and the warm phase of the El Niño are significant at a = .05, but the autoregressive coefficient is not.

Table  shows the parameter estimates for a Poisson GLM with dispersion after modifying the time-series model by dividing the basin into grids, and adding the offset factors for latitude and longitude. The associated t-values on the model coefficients indicate that every factor is significant at a = .05. The fact that the lag-one autoregressive coefficient is significant is unexpected, since it was insignificant in the time-series model. Variances for the global predictors are significantly smaller than the variance of the time series model. Thus, the t-values on the global predictors in the second model are about three times those in the first model.

The conclusions about the parameters and their standard errors are not entirely valid, as the GLM assumes no instantaneous coupling; that is, the conditional distribution given the past is independent for each site. This assumption does not hold because our data consist of hurricanes that pass between adjacent grid boxes creating correlations in annual counts between adjacent boxes. This suggests adding spatial structure to the model in the form of coupling parameters. Note that the GLM formulation does not require temporal independence, since the assumed MLE in the GLM estimator is the same as the actual MLE for an autoregressive time series. Thus, in the absence of instantaneous coupling, we can use the Poisson GLM treating past values of the response as covariates in the model.

Model estimates for the final model are derived using the results of a pseudo likelihood estimator as inputs to the MCML estimator. The MCML estimator is then run iteratively with 1000 samples at each stage. Table  shows the MCML estimates after four iterations of the estimator. Parameter estimates converge with changes between the third and fourth iteration of less than 0.6 times the estimated standard error. The root mean square change is 0.20 s.

Results from the full TPSTAR model indicate all couplings are positive and significant. Note again that the lag-one coupling, which was not significant in the time-series model, is significant in the space-time model. This is a new finding that provides evidence for hurricane path persistence over successive years. Locations that were threatened by a hurricane one year are more likely to be threatened again in the next year. As with the second model, the climate predictors are significant.

Estimates from the MCML estimator appear to be reasonable. The model takes into consideration both spatial and temporal couplings. For example, the parameter estimates are the same sign, but smaller in the final model with instantaneous couplings as compared to the Poisson GLM without instantaneous couplings, or the time-series model. This makes sense, since a positive coupling causes the expected value of any statistic to be more sensitive to changes in the predictor value than would be expected in the absence of the coupling. The parameters and the estimated standard errors from the modified MCML estimator are smaller than those obtained with the first two models. This reduction in standard error might be real owing to the addition of instantaneous coupling and offsets, but it is more likely an artifact of the apparent increase in total hurricane counts from 435 in the time series model to 1676 in the TPSTAR model.

Table  shows the estimated correlations between the coupling parameters and model intercept. The correlations between the temporal and spatial coupling parameters are small, whereas the correlation between the east-west and north-south coupling parameters are considerably larger. Surprisingly the intercept is not strongly correlated to the coupling parameters.

5. Model selection

Backward elimination is applied to the TPSTAR model to arrive at a final model of the spatial/temporal variations in seasonal hurricane activity. The procedure makes use of the estimated changes in the Schwartz's Bayesian Information Criteria (SBC) and the Akaike Information Criteria (AIC). In fact, we need to use changes in criterion levels because the TPSTAR model generates rough estimates of the deviance or predicted values. Thus, common statistics for model selection cannot be calculated accurately. The deviance requires knowledge of the normalizing constant for the distribution, which cannot be estimated accurately in the presence of strong coupling.

We test significance with the Wald test, and use this to compare the difference in SBC or AIC between two models. Both AIC and SBC are used for model selection. Although no asymptotes for either of these statistics exist for our model, for a stationary Gaussian time series, minimizing the AIC gives the model with the smallest predictive error, whereas minimizing the SBC gives a consistent model (Brockwell 1991).

Assume we have two nested models with the second model generated from the first by adding an additional factor with q levels, for a total of p parameters. Let q=[q₁, ¼, q_q] be the parameter vector with q Î R^p and q_q Î R^q for the added factor. Let S_q be the covariance matrix for the added factor, then the Wald test is H₀: q_q=0 and the statistic W = q_q¢S_q^-1q_q, and asymptotically in n, W has a c² distribution with q degrees of freedom. Since AIC=-2logl((q)) + 2p and the asymptotic distributions for W and 2logl(q)-2logl(q₁) are equal, we can approximate the change in AIC for adding the factor with q levels as DAIC_F » -W+2·q. Since SBC=-2logl(q)) +p·log(n), we can calculate the change in SBC for adding a factor with q levels as DSBC_F » -W+q·log(n). Now, if we remove a factor with q levels then the changes in the information are DAIC_B » W-2·q and DSBC_B » W-q·log(n).

Let us consider four versions of the TPSTAR model with instantaneous coupling. The first version is without offsets, so it has three coupling parameters, the intercept, and five climate covariates giving a total of p₁=9 parameters. Version two adds the longitude factor increasing the number of parameters to p₂=p₁+8=17. Version three adds the latitude factor to model one increasing the number of parameters to p₃=p₁+8=14. Version four adds both factors to the original model for a total of p₄=p₁+8+5=22 parameters. First, we test the estimated change in SBC excluding the latitude factor or longitude factor. If the change is negative for either factor we remove this factor from the model. Then we rerun the model and test the change in SBC and AIC by removing the other offset factor. Table  shows these results, which indicate removing the latitude factor based on SBC. Results for the final selected model are shown in Table .

Removal of the latitude factor reduces the global covariate parameters to a small degree, while increasing g_1c. The parameter variances also increase slightly. These results are expected since we are moving the variance explained by the latitude factor into the unexplained variance, which shows as an increase in parameter variance. Note that although removing the latitude factor reduces the SBC it increases the AIC, so we keep the full model when using it to make forecasts.

6. Model fit issues and a hindcast

a. strong coupling

b. grid size

We experimented with larger spatial neighborhoods by adding the current period's four diagonal neighbor sites using a single parameter to the model. These sites are the NE, SE, SW and NW neighbor regions, with parameter g_0d (the diagonal term). We have the same model described by Eq. 2, with Eq. 3 changed to

log(lt,ij) = aij(t)+ g0h(Ht,i+1,j+Ht,i-1,j)+g0v (Ht,i,j+1+Ht,i,j-1)+ g0d (Ht,i+1,j+1+ Ht,i-1,j+1+Ht,i-1,j-1+Ht,i+1,j-1) + g1c (Ht-1,i,j)

and the local neighborhood of site (t,ij), h_¶(t,ij), is the vector of nine response values.

We run the modified MCML estimator for several iterations, using the previous values of the parameter estimates from the full model in Table , and an initial MPLE estimate for g_0d of 0.0 We estimate that g_0d=0.090,s_g0d=0.021 ,t=4.3 . Adding the diagonal parameter to the model significantly affects only the estimates for g_0h and g_0v reducing them to 0.319 and 0.209 respectively, while increasing their standard errors to 0.017 and 0.026, respectively. The correlation matrix of the coupling parameters and the intercept show that while the other terms have correlations less than 0.07 with each other g_0d is significantly correlated. This term is not only significant, but improves the model. If this term is added both the estimated AIC and SBC are reduced by 17.2 and 10.9 respectively. Thus, future models should consider both larger spatial and temporal neighborhoods.

There are concerns about increasing the size of the neighborhood. If the new couplings are positive the estimator may fail to converge. Because we have more site values, the model will required additional boundary values and parameters. One may reduce the number of parameters by using combinations of various canonical parameters. For example, the parameter g_0d is really the sum of two parameters one for the NW-SE diagonal and another for the NE-SW diagonal.

c. hindcast

The model can be used to forecast hurricane activity when lagged values of the covariates are included. At each stage of the MCML estimator, we generate samples from the distribution of X_t conditioned on {X_s  :  s < t } using the observations { X_t-1,¼X_t-p }, the covariates at time t, z_k(t),   k=1¼P the model parameters and estimated boundary values. For our application we used the full model parameters given in Table , so as to have the smallest prediction error. Also, we estimate the boundary values using the mean for each cell.

As a single test case, we generate 103 sample forecasts of hurricane activity for 1994 in each 6^° by 6^° region of grid S using the 1994 values of the five global covariates. Hurricane activity during 1994, which was below the long-term average, is not included at any phase of model development. The spatial distribution and intensity of hindcast values are plotted in Fig. . Table  shows the mean and standard deviation for each region along with a comparison to the actual 1994 values for that region. In this table, if at least one 1994 hurricane entered the grid box the actual value is one. Results indicate reasonable agreement with observations. In general, the hindcast indicates a greater probability in grid boxes that were actually hit.

Results are summarized in Fig.  which shows boxplots of the predicted probability as a function of actual occurrence. Seven grid boxes were affected by hurricanes during 1994. The average hindcast probability over those sites is 0.66 with a median value of 0.59. This compares with an average and median probability of 0.43 and 0.41, respectively for the 33 sites not affected by hurricanes during 1994. Overall this single case study supports the contention that the TPSTAR model, or similar spatial count models, might be useful tools in predicting regional hurricane activity over the North Atlantic basin.

7. Summary and conclusions

We introduce and apply a space-time count process model to North Atlantic hurricane activity. The model uses the best-track data consisting of historical hurricane positions and intensities together with climate variables to determine local space-time coefficients of a truncated Poisson process. The model, referred to as a truncated Poisson space-time autoregressive (or TPSTAR) model, is motivated by first examining a time-series model for the entire domain. Then a Poisson generalized linear model is considered that uses grids boxes within the domain and adds offset factors for latitude and longitude. A natural extension is then made that includes instantaneous local and autoregressive coupling between the grids. A final version of the model is found by backward selection of the predictors based on values of SBC and AIC. A single hindcast is performed on the 1994 hurricane season using a model having five nearest neighbors and statistically significant couplings. The parameters in the TPSTAR model are estimated using MPLE. The model showed promise as a potential forecast tool.

Several conclusions concerning the application of the TPSTAR model to seasonal North Atlantic hurricane activity are reached:

• Seasonal hurricane activity over the North Atlantic basin can be modeled as a space-time Poisson process.
• There appears to be some temporal autocorrelation in hurricane activity between seasons.
• Dividing the region into grid boxes and adding coupling increases the significance of the global predictors.
• Increasing the neighborhood size can improve the model at the cost of estimating additional parameters and boundary values.
• From the MCML estimates t-tests can be performed to determine if the instantaneous couplings are significant. If they are not we can use the simpler Poisson GLM model.
• Grid size is important to the model formulation. The spatial parameters g_0v and g_0h do not scale with changing grid size, and the MCML estimator may fail to converge for small grid size.
• A hindcast case study provides no evidence against the hypothesis that the TPSTAR modeling procedure can be useful as a climate prediction tool.

Acknowledgements. The study was partially funded by the National Science Foundation (ATM-0086958) and the Risk Prediction Initiative of the Bermuda Biological Station for Research (RPI-99-001).

References

• Arguez A, Elsner JB. 2001. Trends in U.S. tropical cyclone mortality during the 20th century. The Florida Geographer. in press.
• Bartlett MS. 1968. A further note on nearest-neighbor models. Journal of the Royal Statistical Society. 131A: 579-580.
• Besag JE. 1974. Spatial interaction and the statistical analysis of lattice systems. Journal of the Royal Statistical Society. 36B: 192-225.
• Brockwell P, Davis A. 1991. Time Series: Theory and Methods, Second Edition. Springer-Verlag: New York.
• Elsner JB, Kara AB. 1999. Hurricanes of the North Atlantic: Climate and Society. Oxford University Press: New York.
• Elsner JB, Bossak BH. 2001. Changes in hurricane landfall probabilities in response to variations in climate. Hurricanes and Typhoons: Past, Present, and Future. Murnane RJ, Liu Kb. Eds. Columbia University Press. in preparation.
• Elsner JB, Kara AB, Owens MA. 1999. Fluctuations in North Atlantic hurricanes. Journal of Climate. 12: 427-437.
• Elsner JB, Lehmiller GS, Kimberlain TB. 1996. Early August forecasts of Atlantic tropical storm activity for the balance of the 1996 season, using Poisson models. Experimental Long-lead Forecast Bulletin. 5(3): 26-28.
• Elsner JB, Niu XF, Tsonis AA. 1998. Multi-year prediction model of North Atlantic hurricane activity. Meteorology and Atmospheric Physics. 68: 42-51.
• Elsner JB, Liu Kb, Kocher B. 2000. Spatial variations in major U.S. hurricane activity: Statistics and a physical mechanism. Journal of Climate. 13: 2293-2305.
• Geyer CM. 1994. On the convergence of Monte Carlo maximum likelihood calculations. Journal of the Royal Statistical Society. 56B: 261-274.
• Gilks WR, Richardson S, Spiegelhalter DJ. 1996. Markov Chain Monte Carlo in Practice. Chapman & Hall: London.
• Gray WM. 1984. Atlantic seasonal hurricane frequency: Part II. Forecasting its variability. Monthly Weather Review. 112: 1669-1683.
• Jagger TH. 2000. ``A Truncated Space-Time Autoregressive Model." Ph.D. Thesis. Department of Statistics, Florida State University.
• Jagger TH, Elsner JB, Niu XF. 2001. A dynamic probability model of hurricane winds in coastal counties of the United States. Journal of Applied Meteorology. 40: 853-863.
• Landsea CW, Gray WM. 1992. The strong association between Western Sahel monsoon rainfall and intense Atlantic hurricanes. Journal of Climate. 5: 435-453.
• Lehmiller GS, Kimberlain TB, Elsner JB. 1997. Seasonal prediction models for North Atlantic basin hurricane location. Monthly Weather Review. 125: 1780-1791.
• Neumann CJ, Jarvinen BR, McAdie CJ, Hammer GR. 1999. Tropical Cyclones of the North Atlantic Ocean, 1871-1998. National Oceanic and Atmospheric Administration.
• Sharpiro LJ. 1987. Month-to-month variability of the Atlantic tropical circulation and its relationship to tropical storm formation. Monthly Weather Review. 115: 2598-2614.
• Whittle P. 1963. Stochastic processes in several dimensions. Bulletin of the International Statistical Institute. 40: Book 2.
• Winkler G. 1995. Image Analysis, Random Fields, and Dynamic Monte Carlo Methods, Springer-Verlag: Berlin.
• Wu H. 1994. ``Regression Models for Spatial Binary Data with Application to the Distribution of Plant Species.'' Ph.D. Thesis. Department of Statistics, Florida State University.