#Landfall hurricanes 1851-2005 (Source = HURDAT)
year<-c(1851:2005)

#1851-2005 data from Willis E @ Climate Audit
storms.all<-c(6,5,8,5,5,6,4,6,8,7,8,6,9,5,7,7,9,4,10,11,8,5,5,7,6,5,8,12,8,11,7,6,4,4,8,12,19,9,9,4,10,9,12,7,6,7,6,11,9,7,12,5,10,5,5,11,5,10,11,5,6,7,6,1,5,14,3,5,3,4,6,4,7,8,2,11,7,6,3,2,9,11,21,11,6,16,9,8,5,8,6,10,10,11,11,6,9,9,13,13,10,7,14,10,13,8,8,10,11,7,11,5,9,12,6,11,8,8,18,10,13,7,8,11,9,10,6,12,9,11,12,6,4,13,11,6,7,12,11,14,8,7,8,7,19,13,8,14,12,15,15,12,16,15,28)
hurricanes.landfall<-c(2,4,1,3,1,2,1,1,1,3,3,0,0,0,2,1,2,0,4,3,3,0,2,1,1,2,2,2,3,4,2,3,1,0,1,7,4,4,1,0,2,0,5,2,1,3,1,3,3,1,2,0,2,2,0,4,0,2,5,2,2,2,3,0,3,6,1,1,1,2,2,0,1,2,1,3,0,2,2,0,0,2,5,2,2,3,0,2,1,2,2,2,1,3,3,1,3,3,3,3,0,1,3,3,3,1,1,0,3,2,1,0,1,4,1,2,1,1,2,1,3,1,0,1,1,1,1,0,3,1,0,0,1,1,6,2,1,1,3,0,1,1,1,0,2,2,1,3,3,0,0,1,2,6,6)

#####################################
# Regression analysis of storm data #
#####################################

#all storms
summary(lm(storms.all~year))                                #all data
1851-2995
summary(allstorms.lm<-lm(storms.all[124:155]~c(1974:2005))) #subset
1974-2005
allstorms.lm.predict<-predict(allstorms.lm,level=.95,interval="confidence",se.fit=T)

#landfall hurricanes
summary(lm(hurricanes.landfall~year))
summary(landfall.lm<-lm(hurricanes.landfall[124:155]~c(1974:2005)))
landfall.lm.predict<-predict(landfall.lm,level=.95,interval="confidence",se.fit=T)

###################
# output analysis #
###################
win.graph(height=10,width=8)

#set graphic parameters

par(mfcol=c(3,2))
par(las=1)
par(cex=0.7)

#landfall hurricanes (in graphic column 1)
plot(year[124:155],storms.all[124:155],type="l",xlab="",ylab="Count (all
storms)")
points(year[124:155],storms.all[124:155],pch=21)
lines(year[124:155],allstorms.lm.predict$fit[,1])
lines(year[124:155],allstorms.lm.predict$fit[,2])
lines(year[124:155],allstorms.lm.predict$fit[,3])
text(1982,25,"(r.sq=0.26, p=0.0015)")
plot(year,storms.all,type="l",xlab="",ylab="Count (all storms)")
pacf(storms.all,main="",lwd=4)

#All storms (in graphic column 2)
plot(year[124:155],hurricanes.landfall[124:155],type="l",xlab="",ylab="landfalling
hurricanes")
points(year[124:155],hurricanes.landfall[124:155],pch=21)
lines(year[124:155],landfall.lm.predict$fit[,1])
lines(year[124:155],landfall.lm.predict$fit[,2])
lines(year[124:155],landfall.lm.predict$fit[,3])
text(1979,5.5,"trend n.s.\n(r.sq=0.08, p=0.08)")
plot(year,hurricanes.landfall,type="l",xlab="",ylab="landfalling
hurricanes")
pacf(hurricanes.landfall,main="",lwd=4)