\chapter{Impact Models}
\label{chap:impactmodels}
\SweaveOpts{keep.source=TRUE, pdf=FALSE, prefix.string=Chap13, grdevice=tikz.Swd}
<<options and clear workspace, echo=FALSE, eval=TRUE>>=
options(digits=3, show.signif.stars=FALSE, width=53)
rm(list=ls())
require(tikzDevice)
source("../SweaveTikZ.R")
@
%data files used: source("Losses.txt"); load("best.use.RData"); load("ncdataframe.RData")
%required packages: maps, ggplot2, oce, sp, quantreg
%source code: source("LossGui.R"); source("plot.windrose2.R"); source("getmax.R"); load("SavedObjects.R")

\begin{quote} ``The skill of writing is to create a context in which other people can think.''
\end{quote}
\indent---Edwin Schlossberg\\

In this chapter we show broader applications of our models and methods.  We focus on impact models.  Hurricanes are capable of generating large financial losses.  We begin with a model that estimates extreme losses conditional on climate covariates.  We then describe a method for quantifying the relative change in potential losses over several decades.

\section{Extreme Losses}

Financial losses are directly related to fluctuations in hurricane climate.  Environmental factors influence the frequency and intensity of hurricanes at the coast (see Chapter~\ref{chap:frequencymodels} and \ref{chap:intensitymodels}).  Thus it is not surprising that these same environmental signals appear in estimates of total losses.

Economic damage is the loss associated with a hurricane's direct impact.\footnote{Direct impact losses do not include losses from business interruption or other macroeconomic effects including demand surge and mitigation.}  A normalization procedure adjusts the loss estimate to what it would be if the hurricane struck in a recent year by accounting for inflation and changes in wealth and population, plus a factor to account for changes in the number of housing units exceeding population growth.  The method produces hurricane loss estimates that can be compared over time \citep{PielkeEtAl2008}.

\subsection{Exploratory analysis}

You focus on losses exceeding one billion (\$ U.S.) that have been adjusted to 2005. The loss data are available in {\it Losses.txt} in JAGS format (see Chapter~\ref{chap:spatialmodels}).  Input the data by typing
<<input data, echo=TRUE, eval=TRUE>>=
source("Losses.txt")
@
The log-transformed loss amounts are in \verb@y@.  The annual number of loss events are in \verb@L@.  The data cover the period 1900--2007.  More details about these data are given in \cite{JaggerEtAl2011}.  You begin by plotting a time series of the number of losses and a histogram of total loss per event.
<<time series and histogram, fig=FALSE, echo=TRUE, eval=FALSE>>=
layout(matrix(c(1, 2), 1, 2, byrow=TRUE),
   widths=c(3/5, 2/5))
plot(1900:2007, L, type="h", xlab="Year",
   ylab="Number of Loss Events")
grid()
mtext("a", side=3, line=1, adj=0, cex=1.1)
hist(y, xlab="Loss Amount ($ log)",
   ylab="Frequency", main="")
mtext("b", side=3, line=1, adj=0, cex=1.1)
@
\begin{figure}
\centering
<<losseventsamount, fig=TRUE, echo=FALSE, eval=TRUE, width=4.6, height=3>>=
par(las=1, mgp=c(2, .4, 0), tcl=-.3)
layout(matrix(c(1, 2), 1, 2, byrow=TRUE), widths=c(.57, .43))
plot(1900:2007, L, type="h", xlab="Year", ylab="Number of loss events")
grid()
points(1900:2007, L, type="h")
mtext("a", side=3, line=1, adj=0, cex=1.1)
hist(y, main="", ylab="Frequency",
  xlab="Loss amount [log$_{10}$ 2005 \\$]")
mtext("b", side=3, line=1, adj=0, cex=1.1)
@
\vspace{-.75cm}
\caption{Loss events and amounts. (a) Number of events and (b) amount of loss.}
\label{fig:losseventsamounts}
\end{figure}
Plots are shown in Fig.~\ref{fig:losseventsamounts}.  The annual number of loss events varies between \Sexpr{min(L)} and \Sexpr{max(L)}.  There appears to be slight increase in the number over time.  The frequency of extreme losses per event decreases nearly linearly on the log scale.  A time series of the amount of loss indicates no trend.  The mean value is the expected loss while the standard deviation is associated with the unexpected loss.  Tail values are used to estimate the `value-at-risk' (VaR) in financial instruments used by the insurance industry.

\subsection{Conditional losses}

You assume a Poisson distribution adequately quantifies the number of extreme loss events\footnote{There is no evidence of loss event clustering.} and a generalized Pareto distribution (GPD) quantifies the amount of losses above the billion dollar threshold (see Chapter~\ref{chap:intensitymodels}).  The frequency of loss events and the amount of losses given an event vary with your covariates including sea-surface temperature (SST), the Southern Oscillation Index (SOI), the North Atlantic Oscillation (NAO) and sunspot number (SSN) as described in Chapter~\ref{chap:datasets}.

The model is written in JAGS code available in {\it JAGSmodel3.txt}.  More details about the model are available in \cite{JaggerEtAl2011}.  Posterior samples from the model (see Chapter~\ref{chap:bayesianmodels}) are available through a graphical user interface (GUI). Start the GUI by typing
<<source the model gui, fig=FALSE, echo=TRUE, eval=FALSE>>=
source("LossGui.R")
@
If you are using a MAC you will need to download and install the GTK+ framework.  You will also probably need the package {\bf gWidgetstcltk}.  Use the slider bars to vary the covariates.  The covariates are scaled between $\pm$2 s.d. Select OK to bring up the return level graph.  You can select the graph to appear as a bar or dot plot.

The GUI allows you to easily ask `what if' questions related to future damage losses.  As one example, Fig.~\ref{fig:posteriorreturnlevels} shows the posterior predictive distributions of return levels for individual loss events using two climate scenarios.  For each return period, the black square shows the median and the red squares show the 0.05, 0.25, 0.75, and 0.9 quantile values from the posterior samples.  The left panel shows the losses when SST = $-$0.243$^\circ$C, NAO = $+$.7 s.d., SOI = $-$1.1 s.d., and SSN = 115. The right panel shows the losses when SST = $+$0.268$^\circ$C, NAO = $-$1.4 s.d., SOI = $-$1.1 s.d., and SSN = 9.
<<generate posterior samples, fig=FALSE, echo=FALSE, eval=TRUE>>=
load("SavedObjects.R")
rp = c(5, 10, 20, 50, 100, 200, 500)
cov1 = data.frame(ssn = 115, soi = -1.1, nao = .7, sst = -.243)
post1 = jagsgpdposterior(newcov=cov1, yp=rp, x=jagSampleMatrix, thin=10,  N=NULL, u=9,
  mask=parmaskgpd, yrange=c(8, 16), yl=18, type="Total Loss", mfrow=c(1, 1),
  etafun=c(exp, exp, identity))
cov2 = data.frame(ssn = 9, soi = -1.1, nao = -1.4, sst = .268)
post2 = jagsgpdposterior(newcov=cov2, yp=rp, x=jagSampleMatrix, thin=10,  N=NULL, u=9,
  mask=parmaskgpd, yrange=c(8, 16), yl=18, type="Total Loss", mfrow=c(1, 1),
  etafun=c(exp, exp, identity))
@
\begin{figure}
\centering
<<posteriorreturnlevels, fig=TRUE, echo=FALSE, eval=TRUE, width=4.6, height=3.3>>=
clrs = c("red", "darkred", "black", "darkred", "red")
qtls = c(3, 4, 5, 6, 7)  # 5, 25, 50, 75, 90 percentiles
par(mfrow=c(1, 2), las=1, mgp=c(2, .4, 0), tcl=-.3)
plot(rp, rep(9, 7), xlim=c(3, 700), ylim=c(9, 14), type="n", bty="n", log="x",
  ylab="Return level [log$_{10}$ 2005 \\$]",
  xlab="Return period [yr]",
  cex.axis=1, cex.lab=1)
grid(equilogs=FALSE)
for(i in 1:length(rp)){
points(rep(rp[i], 5), post1$quantiles[qtls, 1, i], col=clrs, pch=15, cex=0.8)
}
mtext("a", side=3, line=0, adj=0, cex=1.1)
plot(rp, rep(9, 7), xlim=c(3, 700), ylim=c(9, 14), type="n", bty="n", log="x",
  ylab="Return level [log$_{10}$ 2005 \\$]",
  xlab="Return period [yr]",
  cex.axis=1, cex.lab=1)
grid(equilogs=FALSE)
for(i in 1:length(rp)){
points(rep(rp[i], 5), post2$quantiles[qtls, 1, i], col=clrs, pch=15, cex=0.8)
}
mtext("b", side=3, line=0, adj=0, cex=1.1)
@
\vspace{-.75cm}
\caption{Loss return levels.  (a) weak and few versus (b) strong and more.}
\label{fig:posteriorreturnlevels}
\end{figure}
The first scenario is characterized by covariates related to fewer and weaker hurricanes and the second scenario is characterized by covariates related to more and stronger hurricanes.  The loss distribution changes substantially. Under the first scenario, the median return level of a 50-year loss event is \$\Sexpr{round(10^post1$quantiles[5, 1, 4]/10^9,1)} bn; this compares with a median return level of a 50-year loss event of \$\Sexpr{round(10^post2$quantiles[5, 1, 4]/10^9,1)} bn under the second scenario.

The results are interpreted as a return level for a return period of 50 years with the covariate values as extreme or more extreme than the ones given (about one s.d. in each case).  With four independent covariates and an annual probability of about 16\% that a covariate is more than one standard deviation from the mean, the chance that all covariates will be this extreme or more in a given year is less than 0.1\%.  The direction of change is consistent with a change in U.S. hurricane activity using the same scenarios and inline with your understanding of hurricane climate variability in this part of the world.

The model is developed using aggregate loss data for the entire United States susceptible to Atlantic hurricanes.  It is possible to model data representing a subset of losses capturing a particular insurance portfolio, for example.  Moreover since the model uses MCMC sampling it can be extended to include error estimates on the losses.  They can also make use of censored data where you know that losses exceeded a certain level but you do not have information about the actual loss amount.

\subsection{Industry loss models}

Hazard risk affects the profits and losses of companines in the insurance industry. Some of this risk is transferred to the performance of securities traded in financial markets.  This implies that early reliable information from loss models is useful to investors.

A loss model consists of three components; cyclone frequency and intensity, vulnerability of structures, and loss distributions.  They are built using a combination of numerical, empirical, and statistical approaches.  The intensity and frequency components rely on expanding the historical cyclone data set.  This is usually done by resampling and mixing cyclone attributes to generate thousands of synthetic cyclones.

The approach is useful for estimating probable maximum losses (PML), which is the expected value of all losses greater than some high value.  The value is chosen based on the quantile of the loss say $L(\tau)$, where $\tau=1-1/\hbox{RP}$ and where RP is the return period of an event.  The PML can be estimated on a local exposure portfolio or a portfolio spread across portions of the coast.  The PML is used to define the amount of money that a reinsurer or primary insurer should have available to cover a given portfolio. In order to estimate the 1-in-100 year PML a catalogue that contains at least an order of magnitude more synthetic cyclones than contained in the historical data set is required.

Your model above uses only the set of historical losses.  By doing so it allows you to anticipate futures losses on the seasonal and multi-year time scales without relying on a catalogue of synthetic cyclones.  The limitation however is that the losses must be aggregated over a region large enough to capture a sufficient number of loss events.  The aggregation can be geographical, like Florida, or across a large and diverse enough set of exposures.

Your model is also limited by the quality of the historical loss data. Although the per storm damage losses have been adjusted for increases in coastal population the number of loss events have not.  A cyclone making landfall early in the record in a region void of buildings did not generate losses, so there is nothing to adjust.  This is not a problem if losses from historical hurricanes are estimated using constant building exposure data.

\section{Future Wind Damage}

A critical variable in a loss model is hurricane intensity.  Here we show you a way to adjust this variable for climate change.  The methodology focuses on a statistical model for cyclone intensity trends, the output of which can be used as input to a hazard model.  Details are found in \cite{ElsnerEtAl2011}.

\subsection{Historical catalogue}

You begin by determining the historical cyclones relevant to your location of interest.  This is done with the \verb@get.tracks@ function in {\bf getTracks.R} (see Chapter~\ref{chap:datasets}).  Here your interest is a location on Eglin Air Force Base (EAFB) with a latitude 30.4$^\circ$N and longitude $-$86.8$^\circ$E.  You choose a search radius of 100 nmi as a compromise between having enough cyclones to fit a model and having only those that are close enough to cause damage.  You also specify a minimum intensity of 33~m~s$^{-1}$.
<<get locally relevant tracks, echo=TRUE, eval=TRUE>>=
load("best.use.RData")
source("getTracks.R")
lo = -86.8; la = 30.4; r1 = 100
loc = data.frame(lon=lo, lat=la, R=r1)
eafb = get.tracks(x=best.use, locations=loc,
  umin=64, N=200)
@
Tracks meeting the criteria are given in the list object \verb@eafb$tracks@ with each component a data frame containing the attributes of individual cyclones from \verb@best.use@ (see Chapter~\ref{chap:datasets}).

Your interest is restricted further to a segment of each track near the coast.  That is, you subset your tracks keeping only the records when the cyclone is near your location.  This is done using the \verb@selectTrackSubset@ function and specifying a radius for clipping the tracks beyond 300 nmi.  Convert the translation speed (\verb@maguv@) to m~s$^{-1}$.
<<segment the tracks, echo=TRUE, eval=TRUE>>=
r2 = 300
eafb.use = getTrackSubset(tracks=eafb$tracks,
  lon=lo, lat=la, radius=r2)
eafb.use$maguv = eafb.use$maguv * .5144
@
The output is a reduced data frame containing cyclone locations and attributes for track segments corresponding to your location.  Finally you remove tracks that have fewer than 24 hr of attributes.
<<remove short tracks, echo=TRUE, eval=TRUE>>=
x = table(eafb.use$Sid)
keep = as.integer(names(x[x >= 24]))
eafb.use = subset(eafb.use, Sid %in% keep)
@
You plot the tracks on a map (Fig.~\ref{fig:localtracksmap}) reusing your code from Chapter~\ref{chap:datasets}.
\begin{figure}
\centering
<<localtracksmap, fig=TRUE, echo=FALSE, eval=TRUE, width=3.5, height=3.5>>=
t1 = split(eafb.use, eafb.use$Sid)
require(maps)
par(las=1, mgp=c(2, .4, 0), tcl=-.3)
map("world", ylim=c(22, 37), xlim=c(-95, -77))
for(i in 1:length(t1)){
Lo = t1[[i]]$lon
La = t1[[i]]$lat
n = length(Lo)
lines(Lo, La, lwd=2)
arrows(Lo[n - 1], La[n - 1], Lo[n], La[n], lwd=2,
   length=.1)
}
points(lo, la, col="red", pch=19)
box()
@
\vspace{-.75cm}
\caption{Tracks of hurricanes affecting EAFB.}
\label{fig:localtracksmap}
\end{figure}
<<ggplot map of tracks, echo=FALSE, eval=FALSE>>=
require(maps)
require(ggplot2)
states = c("texas", "louisiana", "mississippi", "alabama", "florida", "georgia", "south carolina")
ggplot(eafb.use, aes(lon, lat, group = Sid)) +
  borders("state", states) + theme_bw() +
  ylab(expression(paste("Latitude [",degree,"N]", sep=""))) +
  xlab(expression(paste("Longitude [",degree,"E]", sep=""))) +
  geom_path(color="black")
@
The plot shows a uniform spread of cyclones approaching EAFB from the south with an equal number of cyclones passing to the west as passing to the east.

Your catalogue of historical cyclones affecting EAFB contains \Sexpr{length(unique(eafb.use$Sid))} hurricanes.  You summarize various attributes of these hurricanes with plots and summary statistics.  Your interest is on attributes as the cyclone approaches land so you first subset on the land marker (\verb@M@).
<<subset track attributes for sea only, echo=TRUE, eval=TRUE>>=
sea = subset(eafb.use, !M)
@
A graph showing the distributions of translational speed and approach direction is plotted by typing
<<magnitude and direction of approach, fig=FALSE, echo=TRUE, eval=FALSE>>=
par(mfrow=c(1, 2), mar=c(5, 4, 3, 2) + 1, las=1)
hist(sea$maguv, main="",
  xlab="Translational Speed (m/s)")
require(oce)
u = -sea$maguv * sin(sea$diruv * pi/180)
v = -sea$maguv * cos(sea$diruv * pi/180)
wr = as.windrose(u, v, dtheta=15)
plot.windrose(wr, type="count", cex.lab=.5,
  convention="meteorological")
@
\begin{figure}
\centering
<<translationspeed, fig=TRUE, echo=FALSE, eval=TRUE, width=4.6, height=3.3>>=
require(oce)
source("plot.windrose2.R")
u = -sea$maguv * sin(sea$diruv * pi/180)
v = -sea$maguv * cos(sea$diruv * pi/180)
wr = as.windrose(u, v, dtheta=15)
par(mfrow=c(1, 2), las=1, mgp=c(2, .4, 0), tcl=-.3)
hist(sea$maguv, main="", freq=FALSE, xlab="Translational speed [m s$^{-1}$]")
library(MASS)
fit = fitdistr(sea$maguv,"gamma")
sh = fit$estimate[1]
rt = fit$estimate[2]
curve(dgamma(x, shape=sh, rate=rt), from=0, to=30, col="red", lwd=2, add=TRUE)
mtext("a", side=3, line=0, adj=0, cex=1.2)
plot.windrose2(wr, type="count", convention="meteorological")
mtext("b", side=3, line=-1, adj=0, cex=1)
@
\vspace{-.75cm}
\caption{Translation speed and direction of hurricanes approaching EAFB.}
\label{fig:translationspeed}
\end{figure}
The histograms are shown in Fig.~\ref{fig:translationspeed}.  The smoothed curve on the histogram is a gamma density fit using the \verb@fitdistr@ function ({\bf MASS} package).  The wind rose function is from the {\bf oce} package \citep{Kelly2011}. The median forward speed of approaching cyclones is \Sexpr{round(median(sea$maguv),1)}~m~s$^{-1}$ and the most common approach direction is from the southeast.

The correlation between forward speed and cyclone intensity is \Sexpr{round(cor(sea$maguv,sea$WmaxS),2)}.  For the subset of approaching cyclones that are intensifying this relationship tends to get stronger with increasing intensity. The evidence supports the idea of an intensity limit for hurricanes moving too slow due to the feedback of a cold ocean wake.  This is likely to be even more the case in regions where the ocean mixed layer is shallow.

From a broader perspective, Fig.~\ref{fig:intensityvsforwardspeed} shows the relationship between forward speed and cyclone intensity for all hurricanes over the North Atlantic south of 30$^\circ$N latitude moving slower than 12 m~s$^{-1}$.
\begin{figure}
\centering
<<intensityvsforwardspeed, fig=TRUE, echo=FALSE, eval=TRUE, width=4, height=3>>=
load("best.use.RData")
best.use$WmaxS = best.use$WmaxS * .5144
best.use$maguv = best.use$maguv * .5144
tmp.best = subset(best.use, WmaxS >= 33 & lat <= 30 & maguv < 12)
tmp.split = split(tmp.best, tmp.best$Sid)
tmp.wmax = sapply(tmp.split, function(x) max(x$WmaxS))
tmp.minv = sapply(tmp.split, function(x) mean(x$maguv[max(x$WmaxS) - x$WmaxS < 10]))
tmp.df = data.frame(WmaxS=tmp.wmax, FS=tmp.minv)
require(ggplot2)
p1 = ggplot(tmp.df, aes(FS, WmaxS)) +
  xlab("Translation speed [m s$^{-1}$]") +
  ylab("Lifetime maximum intensity [m s$^{-1}$]") +
  stat_bin2d(bins=20) + geom_smooth(color="red") + theme_bw() +
  opts(axis.title.x=theme_text(vjust=0),
       axis.title.y=theme_text(size=10, angle=90),
       strip.background = theme_blank())
print(p1)
#print(p1 + stat_bin2d(bins=20) + geom_smooth() + theme_bw() +
#  stat_quantile(quantiles=c(.25, .5, .75)))
#summary(lm(WmaxS ~ FS, data=tmp.df))
@
\vspace{0cm}
\caption{Lifetime maximum intensity and translation speed.}
\label{fig:intensityvsforwardspeed}
\end{figure}
The plot shows the lifetime maximum intensity as a function of average translation speed, where the averaging is done when cyclone intensity is within 10 m~s$^{-1}$ of its lifetime maximum.

The two-dimensional plane of the scatter plot is binned into rectangles with the number of points in each bin shown on a color scale.  A local regression line (with standard errors) is added to the plot showing the conditional mean hurricane intensity as a function of forward speed.  The line indicates that, on average, intensity increases with speed, especially for slower moving hurricanes \citep{MeiEtal2012}.  The relationship changes sign for cyclones moving faster than about 8 m~s$^{-1}$.  A linear quantile regression (not shown) indicates the relationship is stronger for quantiles above the median although forward speed explains only a small proportion (about 1\%) of lifetime maximum intensity.

\subsection{Gulf of Mexico hurricanes and SST}

Your historical catalogue of \Sexpr{length(unique(eafb.use$Sid))} cyclones is too small to provide an estimate of changes over time.  Instead you examine the set of hurricanes over the entire Gulf of Mexico.  Changes to hurricanes over this wider region will likely be relevant to changes at EAFB.

Subset the cyclones within a latitude by longitude grid covering the region using the period 1900 through 2009, inclusive.
<<subset cyclones within a latitude longitude grid, echo=TRUE, eval=TRUE>>=
llo = -98; rlo = -80
bla = 19; tla = 32
sy = 1900; ey = 2009
gulf.use = subset(best.use, lon >= llo & lon <= rlo
  & lat >= bla & lat <= tla & Yr >= sy & Yr <= ey)
@
Next find the per cyclone maximum intensity for cyclones while in the Gulf of Mexico.
<<find the maximum intensity for all cyclones in the gulf, echo=TRUE, eval=TRUE>>=
source("getmax.R")
GMI.df = get.max(gulf.use, maxfield="WmaxS")
GMI.df$WmaxS = GMI.df$WmaxS
@

You use the July SST over the Gulf of Mexico as your covariate for modeling the changing intensity of Gulf cyclones.  The gridded SST data are in {\it ncdataframe.RData} where the column names are the year and month concatenated as a character string that includes Y and M.  First create a character vector of the column names.
<<character vector of column names for the ncdataframe, echo=TRUE, eval=TRUE>>=
se = sy:ey
cNam = paste("Y", formatC(se, 1, flag="0"), "M07",
  sep="")
@
Then load the data, create a spatial points data frame and extract the July values for North Atlantic.
<<extract the july sst values from the ncdataframe, echo=TRUE, eval=TRUE>>=
load("ncdataframe.RData")
require(sp)
coordinates(ncdataframe) = c("lon", "lat")
sstJuly = ncdataframe[cNam]
@

Next average the SST values in your Gulf of Mexico grid box.  First create a matrix from your vertex points of your grid box.  Then create a spatial polygons object from the matrix and compute the regional average using the \verb@over@ function.  Next make a data frame from the resulting vector and the structure data set of corresponding years.  Finally merge this data frame with your wind data frame from above.
<<gulf sst average, echo=TRUE, eval=TRUE>>=
bb = c(llo, bla, llo, tla, rlo, tla, rlo, bla,
  llo, bla)
Gulfbb = matrix(bb, ncol=2, byrow=TRUE)
Gulf.sp = SpatialPolygons(list(Polygons(list(
  Polygon(Gulfbb)), ID="Gulfbb")))
SST = over(x=Gulf.sp, y=sstJuly, fn=mean)
SST.df = data.frame(Yr=sy:ey, sst=t(SST)[, 1])
GMI.df = merge(SST.df, GMI.df, by="Yr")
@
The data frame has \Sexpr{dim(GMI.df)[1]} rows one for each hurricane in the Gulf of Mexico region corresponding to the fastest wind speed while in the domain.  The spatial distribution favors locations just off the coast and along the eastern boundary of the domain.

\subsection{Intensity changes with SST}

Theory, models, and data provide support for estimating changes to cyclone intensity over time.  The heat-engine theory argues for an increase in the maximum potential intensity of hurricanes with increases in sea-surface temperature.  Climate model projections indicate an increase in average intensity of about 5--10\% globally by the late 21st century with the frequency of the most intense hurricanes likely increasing even more.  Statistical models using a set of homogeneous tropical cyclone winds show the strongest hurricanes getting stronger with increases as high as 20\% per $^\circ$C for the strongest hurricanes.

Your next step fits a model for hurricane intensity change relevant to your catalogue of cyclones affecting EAFB.  The correlation between per cyclone maximum intensity and the corresponding July SST is a mere \Sexpr{round(cor(GMI.df$sst,GMI.df$WmaxS),2)}, but it increases to \Sexpr{round(cor(GMI.df$sst[GMI.df$WmaxS>64],GMI.df$WmaxS[GMI.df$WmaxS>64]),2)} for the set of cyclones above 64~m~s$^{-1}$.

You use a quantile regression model (see Chapter~\ref{chap:intensitymodels}) to account for the relationship between intensity and SST.  Save the wind speed quantiles and run a quantile regression model of lifetime maximum wind speed (intensity) on SST.  Save the trend coefficients and standard errors from the model for each quantile.
<<quantile regression model, echo=TRUE, eval=TRUE>>=
tau = seq(.05, .95, .05)
qW = quantile(GMI.df$WmaxS, probs = tau)
n = length(qW)
require(quantreg)
model = rq(WmaxS ~ sst, data=GMI.df, tau=tau)
trend = coefficients(model)[2, ]
coef = summary(model, se="iid")
ste = numeric()
for (i in 1:n){
  ste[i] = coef[[i]]$coefficients[2, 2]
}
@
Next use a local polynomial regression to model the trend as a change in intensity per $^\circ$C and plot the results.  The regression fit at intensity $w$ is made using points in the neighborhood of $w$ weighted inversely by their distance to $w$. The neighborhood size is a constant of 75\% of the points.
<<model trend as a percent change in intensity per degree C, fig=FALSE, echo=TRUE, eval=TRUE>>=
trend = trend/qW * 100
ste = ste/qW * 100
model2 = loess(trend ~ qW)
pp = predict(model2, se=TRUE)
xx = c(qW, rev(qW))
yy = c(pp$fit + 2 * pp$se.fit,
  rev(c(pp$fit - 2 * pp$se.fit)))
plot(qW, trend, pch=20, ylim=c(-20, 30),
  xlab="Intensity (m/s)",
  ylab="Percent Change (per C)")
polygon(xx, yy, col="gray", border="gray")
for(i in 1:n) segments(qW[i], trend[i] - ste[i],
  qW[i], trend[i] + ste[i])
points(qW, trend, pch=20)
lines(qW, fitted(model2), col="red")
abline(h=0, lty=2)
@
\begin{figure}
\centering
<<trendperdegreec, fig=TRUE, echo=FALSE, eval=TRUE, width=4.6, height=3.3>>=
par(las=1, mgp=c(2, .4, 0), tcl=-.3)
plot(qW, trend, pch=20, ylim=c(-20, 30),
  xlab="Intensity [m s$^{-1}$]",
  ylab="Percent change [/$^\\circ$C]")
polygon(xx, yy, col="gray", border="gray")
for(i in 1:n) segments(qW[i], trend[i] - ste[i], qW[i], trend[i] + ste[i])
points(qW, trend, pch=20)
lines(qW, fitted(model2), col="red")
abline(h=0, lty=2)
@
\vspace{-.75cm}
\caption{Intensity change as a function of SST for Gulf of Mexico hurricanes.}
\label{fig:trendperdegreec}
\end{figure}
Results are shown in Fig.~\ref{fig:trendperdegreec}.    Points are quantile regression coefficients of per cyclone maximum intensity on SST and the vertical bars are the standard errors.  The red line is a local polynomial fit through the points and the gray band is the 95\% confidence band around the predicted fit.  There is little change in intensity for the weaker cyclones but there is a large, and for some quantiles, statistically significant upward trend in intensity for the stronger hurricanes.

\subsection{Stronger hurricanes}

Next you quantify the trend in SST over time using linear regression of SST on year.
<<trend in sst, echo=TRUE, eval=TRUE>>=
model3 = lm(sst ~ Yr, data=SST.df)
model3$coef[2] * 100
@
The upward trend is \Sexpr{round(model3$coef[2]*100,2)}$^\circ$C per century, explaining \Sexpr{round(100*summary(model3)$r.squared)}\% of the variation in July SST over the period of record.  The magnitude of warming you find in the Gulf of Mexico is consistent with reports of between 0.4 and 1$^\circ$C per century for warming of the tropical oceans \citep{DeserEtAl2010}.

Your estimate of the per $^\circ$C SST increase in hurricane intensities (as a function quantile intensity) together with your estimate of the SST warming is used to estimate the increase in wind speeds for each hurricane in your catalogue. You assume your catalogue is a representative sample of the frequency and intensity of future hurricanes, but that the strongest hurricanes will be stronger due to the additional warmth expected. The approach is similar to that used in \cite{MousaviEtAl2011} to estimate the potential impact of hurricane intensification and sea-level rise on coastal flooding.

The equation for increasing wind speeds $w$ of hurricanes in your catalogue of hurricanes affecting EAFB is given by
\begin{equation}
w_{2110} = [1 + \Delta w(q) \cdot \Delta\hbox{SST}] \cdot w
\end{equation}
where $w_{2110}$ is the wind speed one hundred years from now, $\Delta w(q)$ is the fractional change in wind speed per degree change in SST as a function of the quantile speed (red curve), and $\Delta$SST is per century trend in SST.  You certainly do not expect an extrapolation (linear or otherwise) to accurately represent the future but the method provides an estimate of what Gulf of Mexico hurricanes might look like, on average, during the 22nd century.

Additional cyclone vitals including radius to maximum winds, a wind decay parameter, and minimum central pressure need to be added to the catalogue of cyclones to make them useful for storm surge and wind-field components \citep{VickeryEtAl2006} like those included in the U.S. Federal Emergency Management Agency (FEMA) HAZUS model.  Lacking evidence indicating these vitals will change in the future you use historical values.  Otherwise you adopt a method similar to what you used above for cyclone intensity.  In the end you have two cyclone catalogues, one representing the contemporary view and the other representing a view of the future; a view that is consistent with the current evidence and theory of cyclone intensity and which aligns with the consensus view on anthropogenic global warming.

This chapter demonstrated a few ways in which the models and methods described in the earlier are used to answer questions about possible future impacts.  In particular, we showed how data on past insured losses are used to hedge against future losses conditional on the state of the climate.  We also showed how to model potential changes to local hurricane impacts caused by a changing cyclone climatology.

This concludes your study of hurricane climatology using modern statistical methods.  We trust your skills will help you make new reproduceable discoveries in the fascinating field of hurricane climate.  Enjoy.

%END OF BOOK

%NEW STUFF
<<5-year model for proportions using beta regression, echo=FALSE, eval=FALSE>>=
#Medium term rates
load("annual.RData")
annual5year = annual
annualcount = annual[2:31]
annual5yearcount = data.frame(lapply(annualcount,
  function(w) rev(filter(rev(w), c(1, 1, 1, 1, 1),
  sides=1))))
annual5year[2:31] = annual5yearcount

#The sums are y[t]=y[t]+y[t+1]+y[t+2]+y[t+3]+y[t+4]
#The covariates are x[t]

#Sample:
annual5year[1:2, c("Year", "B.1", "B.2", "US.1",
                   "ssn.May")]

#stepwise linear regression on proportion
fullmodel = lm(I(US.1/B.1) ~ soi + sst + ssn + nao, data=annual5year[16:156, ])
step(fullmodel)

model = lm(I(US.1/B.1) ~ sst + ssn, data=annual5year[16:156, ])
summary(model)

model = lm(I(US.1/B.1) ~ sst.Jun + ssn.Jun, data=annual5year[16:156, ])
summary(model)

require(betareg)
logit = betareg(I(US.1/B.1) ~ sst.Jun + ssn.Jun, data=annual5year[16:156, ])

ssnp = seq(0, 201, 1)
sstp = rep(0, length(ssnp))

pre = predict(logit, data.frame(sst.Jun=sstp, ssn.Jun=ssnp))
plot(ssnp, pre, type="l", ylim=c(.1, .8))
points(annual5year$ssn.Jun, annual5year$US.1/annual5year$B.1)

require(ggplot2)
require(grid)
a5 = annual5year[16:156, ]
a5$lp = a5$US.1/a5$B.1
p1 = ggplot(data=a5, aes(sst.Jun, lp)) +  geom_point() + theme_bw()
p1 = p1 + stat_smooth() + ylim(0, .8)
p1 = p1 + xlab(expression(paste("June SST anomaly [", degree, "C]"))) +
  ylab("U.S. landfall proportion over the next 5 years")

p2 = ggplot(data=a5, aes(ssn.Jun, lp)) +  geom_point() + theme_bw()
p2 = p2 + stat_smooth() + ylim(0, .8)
p2 = p2 + xlab("June SSN [count]") + ylab("U.S. landfall proportion over the next 5 years")

grid.newpage()
pushViewport(viewport(layout = grid.layout(1, 2)))
vplayout <- function(x, y)
    viewport(layout.pos.row = x, layout.pos.col = y)
print(p1 +
  opts(title="a", plot.title=theme_text(hjust=0, size=11), aspect.ratio=1),
  vp=vplayout(1, 1))
print(p2 +
  opts(title="b", plot.title=theme_text(hjust=0, size=11), aspect.ratio=1),
  vp=vplayout(1, 2))
@ 

<<LMI versus time relative to local sunset>>=
LMI = read.table("LMI.txt", header=TRUE)

#Sunset time
loc = matrix(c(LMI$lon, LMI$lat), ncol=2)
day = ISOdate(year=LMI$Yr, month=LMI$Mo, day=LMI$Da, tz="GMT")
require(maptools)
sset = sunriset(crds=loc, day, proj4string=CRS("+proj=longlat +datum=WGS84"), 
                direction="sunrise", POSIXct.out=TRUE)

#LMI time
lmit = ISOdatetime(year=LMI$Yr, month=LMI$Mo, day=LMI$Da, hour=LMI$hr, min=0, sec=0, tz="GMT")

#time difference in seconds relative to sunset with negative values before sunset
timediff = difftime(lmit, sset$time, tz="GMT", units="hours")  
thss = as.numeric(timediff, "hours")
thss.df = data.frame(hours=thss)

p0 = ggplot(thss.df, aes(hours)) + geom_histogram(binwidth=1) + theme_bw()
p0 + xlab("Hours relative to local sunrise") + ylab("Frequency of lifetime maximum")

ggplot(thss.df, aes(hours)) + 
   geom_histogram(aes(y = ..count..), colour = 'orange', binwidth = 1) + 
   coord_polar(theta = "x", start = 2*pi, direction = 1, expand = FALSE) 

@