#Set the current working directory to the CA data directory
cwd<-"http://www.climateaudit.org/data/hurricane/"

##################################################
# R-script for Hurricane time-series analysis    #
# by bender, Climate Audit                       #
# Oct 10, 2006                                   #
##################################################

##################################################
# Program works from two datafiles:              #
# 1. Hurdat 1851-2005                            #
# 2. Emanuel PDI from Emanuel (2005)             #
#    (PDI=power dissipation index)               #
##################################################

##################
# Part I. HURDAT #
##################

#All storms
infile1<-paste(cwd,"allstorms.txt",sep="")

allstorms<-read.table(infile1,header=T)
attach(allstorms)

#linear time-trend regression model
summary(lm(count~c(1851:2005)))
summary(allstorms.lm<-lm(count[124:155]~c(1974:2005)))
allstorms.lm.predict<-predict(allstorms.lm,level=.95,interval="confidence",se.fit=T)

#ARMA model
allstorms.arima<-arima(count,order=c(1,0,1))
predict(allstorms.arima,n.ahead=3)
allstorms.arima.5<-arima(count,order=c(5,0,1))
predict(allstorms.arima.5,n.ahead=3)

#Test the effect of an inactive 2006 season
count.06<-c(count,6)
allstorms.arima<-arima(count.06,order=c(5,0,1))
predict(allstorms.arima,n.ahead=3)
summary(allstorms.lm2<-lm(c(count[124:155],6)~c(1974:2006)))
allstorms.lm.2.predict<-predict(allstorms.lm2,level=.95,interval="confidence",se.fit=T)

#Test the effect of an inactive 2007 season
count.07<-c(count,6,6)
allstorms.arima<-arima(count.07,order=c(5,0,1))
predict(allstorms.arima,n.ahead=3)
summary(allstorms.lm2<-lm(c(count[124:155],6,6)~c(1974:2007)))

win.graph(height=10,width=10)
par(mfcol=c(3,2))
par(las=1)
par(cex=0.9)

plot(year[124:155],count[124:155],type="l",xlab="",ylab="Count (all storms)")
points(year[124:155],count[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,count,type="l",xlab="",ylab="Count (all storms)")
lines(year,lowess(count,f=0.2)$y,col="red")
pacf(count,main="",lwd=4)
detach(allstorms)

#Landfalling hurricanes
infile2<-paste(cwd,"landfall.txt",sep="")

landfall<-read.table(infile2,header=T)
attach(landfall)
summary(lm(count~c(1851:2005)))
summary(landfall.lm<-lm(count[124:155]~c(1974:2005)))
landfall.lm.predict<-predict(landfall.lm,level=.95,interval="confidence",se.fit=T)

plot(year[124:155],count[124:155],type="l",xlab="",ylab="landfalling hurricanes")
points(year[124:155],count[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,count,type="l",xlab="",ylab="landfalling hurricanes")
lines(year,lowess(count,f=0.2)$y,col="red")
pacf(count,main="",lwd=4)

detach(landfall)

cor(landfall$count,allstorms$count)
cor(landfall$count[1:123],allstorms$count[1:123])
cor(landfall$count[124:155],allstorms$count[124:155])

cor(landfall$count[1:78],allstorms$count[1:78])
cor(landfall$count[78:155],allstorms$count[78:155])

#time-series comparison (correlation vs. time period)
attach(allstorms)
win.graph(height=6,width=8)
par(las=1)
plot(year,allstorms$count,type="l",ylab="storm count",xlab="",bg="white")
lines(year,landfall$count,col="red")
abline(v=c(1855,1880,1900,1920,1940,1960,1980,2000))
text(1890,27,"1850-1927\nr=0.62")
text(1970,27,"1927-2005\nr=0.49")
legend(1905,25,lty=c(1,1),col=c("black","red"),c("all storms","landfalling hurricanes"),bg="white")

detach(allstorms)

#ARMA(1,1) comparison
arima(allstorms$count,order=c(1,0,1))
arima(landfall$count,order=c(1,0,1))

##################################################
# Part II. Emanuel's (2005) PDI                  #
#                                                #
# Willis Essenbach's new Emanuel (2005) data set #
##################################################

#Emanuel's "1,2,1" filter
filter.em<-c(0.25,0.5,0.25)

infile5<-paste(cwd,"emanuelPDI.txt",sep="")
em<-read.table(infile5,header=T)
attach(em)

win.graph(height=6,width=4)
par(mfcol=c(2,1))
par(las=1)
par(cex=0.8)
hist(AdjPDI)
hist(sqrt(AdjPDI))

summary(AdjPDI.lm<-lm(AdjPDI~Year))
PDI.lm.predict<-predict(AdjPDI.lm,level=.95,interval="confidence",se.fit=T)

#PDI time-trend analysis plot
win.graph(height=4,width=4)
par(las=1)
plot(Year,AdjPDI,type="l",xlab="",ylab="Atlantic PDI")
points(Year,AdjPDI,pch=21)
lines(Year,PDI.lm.predict$fit[,1])
lines(Year,PDI.lm.predict$fit[,2])
lines(Year,PDI.lm.predict$fit[,3])
abline(v=c(1950,1960,1970,1980,1990,2000),col="grey")
text(1977,28,"(r.sq=0.016, p=0.3430)")

#PDI vs SST time-series plot
win.graph(height=6,width=5)
par(mfcol=c(2,1))
par(las=1)
par(cex=0.8)
plot(Year,AdjPDI,type="l",xlab="")
abline(v=c(1950,1960,1970,1980,1990,2000))
plot(Year,HadISST,type="l",xlab="")
abline(v=c(1950,1960,1970,1980,1990,2000))

#SST time-series analysis (for 0x,1x,2x smoothing)
win.graph(height=8,width=9)
par(mfcol=c(3,3))
par(las=1)
par(cex=0.8)
plot(Year,HadISST,type="l",xlab="")
abline(v=c(1950,1960,1970,1980,1990,2000))
pacf(na.omit(HadISST),lwd=5,main="0xS SST")
spectrum(na.omit(HadISST),spans=c(3),main="0xS SST") 

#Do the 1x filtering
HadSST.flt<-filter(HadISST,filter.em)

plot(Year,HadSST.flt,type="l",xlab="")
abline(v=c(1950,1960,1970,1980,1990,2000))
pacf(na.omit(HadSST.flt),lwd=5,main="1xS SST")
spectrum(na.omit(HadSST.flt),spans=c(3),main="1xS SST")

#Do the 2x filtering
HadSST.2.flt<-filter(HadSST.flt,filter.em)

plot(Year,HadSST.2.flt,type="l",xlab="")
abline(v=c(1950,1960,1970,1980,1990,2000))
pacf(na.omit(HadSST.2.flt),lwd=5,main="2xS SST")
spectrum(na.omit(HadSST.2.flt),spans=c(3),main="2xS SST") 

#PDI time-series analysis with 0x,1x,2x smoothing
win.graph(height=8,width=9)
par(mfcol=c(3,3))
par(las=1)
par(cex=0.8)
plot(Year,AdjPDI,type="l",xlab="")
abline(v=c(1950,1960,1970,1980,1990,2000))
#acf(na.omit(AdjPDI),lwd=5)
pacf(na.omit(AdjPDI),lwd=5,main="0xS PDI")
spectrum(na.omit(AdjPDI),spans=c(3),main="0xS PDI")

#Do the 1x filtering
AdjPDI.flt<-filter(AdjPDI,filter.em)

plot(Year,AdjPDI.flt,type="l",xlab="")
abline(v=c(1950,1960,1970,1980,1990,2000))
#acf(na.omit(AdjPDI.flt),lwd=5)
pacf(na.omit(AdjPDI.flt),lwd=5,main="1xS PDI")
spectrum(na.omit(AdjPDI.flt),spans=c(3),main="1xS PDI")

#Do the 2x filtering
AdjPDI.2.flt<-filter(AdjPDI.flt,filter.em)

plot(Year,AdjPDI.2.flt,type="l",xlab="")
abline(v=c(1950,1960,1970,1980,1990,2000))
#acf(na.omit(AdjPDI.2.flt),lwd=5)
pacf(na.omit(AdjPDI.2.flt),lwd=5,main="2xS PDI")
spectrum(na.omit(AdjPDI.2.flt),spans=c(3),main="2xS PDI")

#PDI vs. SST correlation analysis (for 0x,1x,2x smoothing)
cor(HadISST,AdjPDI,use="pairwise.complete.obs")
cor(HadSST.flt,AdjPDI.flt,use="pairwise.complete.obs")
cor(HadSST.2.flt,AdjPDI.2.flt,use="pairwise.complete.obs")

#PDI vs. SST regression analysis (for 0x,1x,2x smoothing)
summary(lm(AdjPDI~HadISST))
summary(lm(AdjPDI.flt~HadSST.flt))
summary(lm(AdjPDI.2.flt~HadSST.2.flt))

#Test various ARMA fits to PDI
arima(AdjPDI,order=c(1,0,0)) #AR(1)
arima(AdjPDI,order=c(0,0,1)) #MA(1)
arima(AdjPDI,order=c(1,0,1)) #ARMA(1,1)
arima(AdjPDI,order=c(1,1,1)) #ARIMA(1,1,1)

####################################################
#time-trend analysis (Mann-Kendall non-parametric) #
####################################################
library(Kendall)

#1974+
MannKendall(AdjPDI[26:58])
MannKendall(AdjPDI.flt[26:58])
MannKendall(AdjPDI.2.flt[26:58])

#1949+
MannKendall(AdjPDI[1:58])
MannKendall(AdjPDI.flt[1:58])
MannKendall(AdjPDI.2.flt[1:58])

################################
# Attribution: Trend vs. cycle #
################################

#spectral coherence
win.graph(height=7,width=5)
par(mfcol=c(3,1))
x<-cbind(as.ts(HadISST[1:55]),as.ts(AdjPDI[1:55]))
x.spc<-spectrum(x,spans=c(3),plot=F)
plot(x.spc,plot.type="coh",main="Emanuel (2005) SST vs. PDI")
x<-cbind(as.ts(HadSST.flt[2:54]),as.ts(AdjPDI.flt[2:54]))
x.spc<-spectrum(x,spans=c(3),plot=F)
plot(x.spc,plot.type="coh",main="Emanuel (2005) SST vs. PDI, 1x smooth")
x<-cbind(as.ts(HadSST.2.flt[3:53]),as.ts(AdjPDI.2.flt[3:53]))
x.spc<-spectrum(x,spans=c(3),plot=F)
plot(x.spc,plot.type="coh",main="Emanuel (2005) SST vs. PDI, 2x smooth")

###################################
# Breakpoint sensitivity analysis #
###################################

#0x smoothing
summary(lm(AdjPDI~HadISST))                             #1949
summary(lm(AdjPDI[11:52]~HadISST[11:52]))               #1959
summary(lm(AdjPDI[21:52]~HadISST[21:52]))               #1969

#1x smoothing
summary(lm(AdjPDI.flt~HadSST.flt))                      #1949
summary(lm(AdjPDI.flt[11:52]~HadSST.flt[11:52]))        #1959
summary(lm(AdjPDI.flt[21:52]~HadSST.flt[21:52]))        #1969

#2x smoothing
summary(AdjPDI.lm<-lm(AdjPDI.2.flt~HadSST.2.flt))       #1949
summary(lm(AdjPDI.2.flt[11:52]~HadSST.2.flt[11:52]))    #1959
summary(lm(AdjPDI.2.flt[21:52]~HadSST.2.flt[21:52]))    #1969

#spectrum of residuals for 2x smoothing
win.graph(height=6,width=4)
par(mfcol=c(2,1))
pacf(residuals(AdjPDI.lm),lwd=5)
spectrum(residuals(AdjPDI.lm),spans=c(3))

#############################################
# What if filtering were done improperly    #
# by preserving endpoints, rather           #
# than eliminating them?                    #
# What differnce do the low 2006 data make? #
#############################################

#################################################################################
# Working with 1949+ data
#################################################################################

#With 2006 data included

#Filter SSTs

#Do the 1x filtering
HadSST.flt<-filter(HadISST,filter.em)                   #filter
HadSST.flt[1]<-HadISST[1]                               #insert lopped-off start point
HadSST.flt[55]<-HadISST[55]                             #insert lopped-off end point (2003)

#Do the 2x filtering
HadSST.2.flt<-filter(HadSST.flt,filter.em)              #filter
HadSST.2.flt[1]<-HadSST.flt[1]                          #insert lopped-off start point
HadSST.2.flt[55]<-HadSST.flt[55]                        #insert lopped-off end point (2003)

#Filter PDIs

#Do the 1x filtering
AdjPDI.flt<-filter(AdjPDI,filter.em)                    #filter
AdjPDI.flt[1]<-AdjPDI[1]                                #insert lopped-off start point
AdjPDI.flt[58]<-AdjPDI[58]                              #insert lopped-off end point (2006)

#Do the 2x filtering
AdjPDI.2.flt<-filter(AdjPDI.flt,filter.em)              #filter
AdjPDI.2.flt[1]<-AdjPDI.flt[1]                          #insert lopped-off start point
AdjPDI.2.flt[58]<-AdjPDI.flt[58]                        #insert lopped-off end point (2006)

summary(lm(AdjPDI~HadISST))
summary(lm(AdjPDI.flt~HadSST.flt))
summary(lm(AdjPDI.2.flt~HadSST.2.flt))

#time-trend analysis (parametric)
summary(AdjPDI.lm<-lm(AdjPDI~Year))
summary(AdjPDI.lm.flt<-lm(AdjPDI.flt~Year))
summary(AdjPDI.lm.2.flt<-lm(AdjPDI.2.flt~Year))

#calculate confidence intervals
PDI.lm.predict<-predict(AdjPDI.lm,level=.95,interval="confidence",se.fit=T)
PDI.flt.lm.predict<-predict(AdjPDI.lm.flt,level=.95,interval="confidence",se.fit=T)
PDI.2.flt.lm.predict<-predict(AdjPDI.lm.2.flt,level=.95,interval="confidence",se.fit=T)

#PDI time-trend analysis plot
win.graph(height=8,width=5)
par(las=1)
par(mfcol=c(3,1))

plot(Year,AdjPDI,type="l",xlab="",ylab="Atlantic PDI, 0x",ylim=c(0,35))
points(Year,AdjPDI,pch=21)
lines(Year,PDI.lm.predict$fit[,1])
lines(Year,PDI.lm.predict$fit[,2])
lines(Year,PDI.lm.predict$fit[,3])
abline(v=c(1950,1960,1970,1980,1990,2000),col="grey")
text(1977,28,"(r.sq=0.016, p=0.3430)")

plot(Year,AdjPDI.flt,type="l",xlab="",ylab="Atlantic PDI, 1x",ylim=c(0,35))
points(Year,AdjPDI.flt,pch=21)
lines(Year,PDI.flt.lm.predict$fit[,1])
lines(Year,PDI.flt.lm.predict$fit[,2])
lines(Year,PDI.flt.lm.predict$fit[,3])
abline(v=c(1950,1960,1970,1980,1990,2000),col="grey")
text(1977,28,"(r.sq=0.030, p=0.1901)")

plot(Year,AdjPDI.2.flt,type="l",xlab="",ylab="Atlantic PDI, 2x",ylim=c(0,35))
points(Year,AdjPDI.2.flt,pch=21)
lines(Year,PDI.2.flt.lm.predict$fit[,1])
lines(Year,PDI.2.flt.lm.predict$fit[,2])
lines(Year,PDI.2.flt.lm.predict$fit[,3])
abline(v=c(1950,1960,1970,1980,1990,2000),col="grey")
text(1977,28,"(r.sq=0.037, p=0.1458)")

#With 2005 data only (2006 PDI excluded)

#Filter PDIs

#Do the 1x filtering
AdjPDI.flt<-filter(AdjPDI[1:57],filter.em)              #filter
AdjPDI.flt[1]<-AdjPDI[1]                                #insert lopped-off start point
AdjPDI.flt[57]<-AdjPDI[57]                              #insert lopped-off end point (2005)

#Do the 2x filtering
AdjPDI.2.flt<-filter(AdjPDI.flt,filter.em)              #filter
AdjPDI.2.flt[1]<-AdjPDI.flt[1]                          #insert lopped-off start point
AdjPDI.2.flt[57]<-AdjPDI.flt[57]                        #insert lopped-off end point (2005)

#time-trend analysis (parametric)
summary(AdjPDI.lm<-lm(AdjPDI[1:57]~Year[1:57]))
summary(AdjPDI.lm.flt<-lm(AdjPDI.flt[1:57]~Year[1:57]))
summary(AdjPDI.lm.2.flt<-lm(AdjPDI.2.flt[1:57]~Year[1:57]))

#calculate confidence intervals
PDI.lm.predict<-predict(AdjPDI.lm,level=.95,interval="confidence",se.fit=T)
PDI.flt.lm.predict<-predict(AdjPDI.lm.flt,level=.95,interval="confidence",se.fit=T)
PDI.2.flt.lm.predict<-predict(AdjPDI.lm.2.flt,level=.95,interval="confidence",se.fit=T)


#PDI time-trend analysis plot
win.graph(height=8,width=5)
par(las=1)
par(mfcol=c(3,1))

plot(Year[1:57],AdjPDI[1:57],type="l",xlab="",ylab="Atlantic PDI, 0x",ylim=c(0,35))
points(Year[1:57],AdjPDI[1:57],pch=21)
lines(Year[1:57],PDI.lm.predict$fit[,1])
lines(Year[1:57],PDI.lm.predict$fit[,2])
lines(Year[1:57],PDI.lm.predict$fit[,3])
abline(v=c(1950,1960,1970,1980,1990,2000),col="grey")
text(1977,28,"(r.sq=0.021, p=0.2847)")

plot(Year[1:57],AdjPDI.flt[1:57],type="l",xlab="",ylab="Atlantic PDI, 1x",ylim=c(0,35))
points(Year[1:57],AdjPDI.flt[1:57],pch=21)
lines(Year[1:57],PDI.flt.lm.predict$fit[,1])
lines(Year[1:57],PDI.flt.lm.predict$fit[,2])
lines(Year[1:57],PDI.flt.lm.predict$fit[,3])
abline(v=c(1950,1960,1970,1980,1990,2000),col="grey")
text(1977,28,"(r.sq=0.048, p=0.1025)")

plot(Year[1:57],AdjPDI.2.flt[1:57],type="l",xlab="",ylab="Atlantic PDI, 2x",ylim=c(0,35))
points(Year[1:57],AdjPDI.2.flt[1:57],pch=21)
lines(Year[1:57],PDI.2.flt.lm.predict$fit[,1])
lines(Year[1:57],PDI.2.flt.lm.predict$fit[,2])
lines(Year[1:57],PDI.2.flt.lm.predict$fit[,3])
abline(v=c(1950,1960,1970,1980,1990,2000),col="grey")
text(1977,28,"(r.sq=0.066, p=0.0491)")

#################################################################################
# As above but with 1974+ data (to Match Emanuel's time-frame) N.B. 1974=Year[26]
#################################################################################

#With 2006 data included

#Filter PDIs

#Do the 1x filtering
AdjPDI.flt<-filter(AdjPDI,filter.em)                    #filter
AdjPDI.flt[26]<-AdjPDI[26]                              #insert lopped-off start point
AdjPDI.flt[58]<-AdjPDI[58]                              #insert lopped-off end point (2006)

#Do the 2x filtering
AdjPDI.2.flt<-filter(AdjPDI.flt,filter.em)              #filter
AdjPDI.2.flt[26]<-AdjPDI.flt[26]                        #insert lopped-off start point
AdjPDI.2.flt[58]<-AdjPDI.flt[58]                        #insert lopped-off end point (2006)

#time-trend analysis (parametric)
summary(AdjPDI.lm<-lm(AdjPDI[26:58]~Year[26:58]))
summary(AdjPDI.lm.flt<-lm(AdjPDI.flt[26:58]~Year[26:58]))
summary(AdjPDI.lm.2.flt<-lm(AdjPDI.2.flt[26:58]~Year[26:58]))

#calculate confidence itnervals
PDI.lm.predict<-predict(AdjPDI.lm,level=.95,interval="confidence",se.fit=T)
PDI.flt.lm.predict<-predict(AdjPDI.lm.flt,level=.95,interval="confidence",se.fit=T)
PDI.2.flt.lm.predict<-predict(AdjPDI.lm.2.flt,level=.95,interval="confidence",se.fit=T)

#PDI time-trend analysis plot
win.graph(height=8,width=5)
par(las=1)
par(mfcol=c(3,1))

plot(Year[26:58],AdjPDI[26:58],type="l",xlab="",ylab="Atlantic PDI, 0x",ylim=c(0,35))
points(Year[26:58],AdjPDI[26:58],pch=21)
lines(Year[26:58],PDI.lm.predict$fit[,1])
lines(Year[26:58],PDI.lm.predict$fit[,2])
lines(Year[26:58],PDI.lm.predict$fit[,3])
abline(v=c(1950,1960,1970,1980,1990,2000),col="grey")
text(1981,28,"(r.sq=0.23, p=0.0045)")

plot(Year[26:58],AdjPDI.flt[26:58],type="l",xlab="",ylab="Atlantic PDI, 1x",ylim=c(0,35))
points(Year[26:58],AdjPDI.flt[26:58],pch=21)
lines(Year[26:58],PDI.flt.lm.predict$fit[,1])
lines(Year[26:58],PDI.flt.lm.predict$fit[,2])
lines(Year[26:58],PDI.flt.lm.predict$fit[,3])
abline(v=c(1950,1960,1970,1980,1990,2000),col="grey")
text(1981,28,"(r.sq=0.35, p=0.0002)")

plot(Year[26:58],AdjPDI.2.flt[26:58],type="l",xlab="",ylab="Atlantic PDI, 2x",ylim=c(0,35))
points(Year[26:58],AdjPDI.2.flt[26:58],pch=21)
lines(Year[26:58],PDI.2.flt.lm.predict$fit[,1])
lines(Year[26:58],PDI.2.flt.lm.predict$fit[,2])
lines(Year[26:58],PDI.2.flt.lm.predict$fit[,3])
abline(v=c(1950,1960,1970,1980,1990,2000),col="grey")
text(1981,28,"(r.sq=0.43, p=0.0000)")

#With 2005 data only (2006 PDI excluded)

#Filter PDIs

#Do the 1x filtering
AdjPDI.flt<-filter(AdjPDI[1:57],filter.em)              #filter
AdjPDI.flt[26]<-AdjPDI[26]                              #insert lopped-off start point
AdjPDI.flt[57]<-AdjPDI[57]                              #insert lopped-off end point (2005)

#Do the 2x filtering
AdjPDI.2.flt<-filter(AdjPDI.flt,filter.em)              #filter
AdjPDI.2.flt[26]<-AdjPDI.flt[26]                        #insert lopped-off start point
AdjPDI.2.flt[57]<-AdjPDI.flt[57]                        #insert lopped-off end point (2005)

#time-trend analysis (parametric)
summary(AdjPDI.lm<-lm(AdjPDI[26:57]~Year[26:57]))
summary(AdjPDI.lm.flt<-lm(AdjPDI.flt[26:57]~Year[26:57]))
summary(AdjPDI.lm.2.flt<-lm(AdjPDI.2.flt[26:57]~Year[26:57]))

#calculate confidence itnervals
PDI.lm.predict<-predict(AdjPDI.lm,level=.95,interval="confidence",se.fit=T)
PDI.flt.lm.predict<-predict(AdjPDI.lm.flt,level=.95,interval="confidence",se.fit=T)
PDI.2.flt.lm.predict<-predict(AdjPDI.lm.2.flt,level=.95,interval="confidence",se.fit=T)

#PDI time-trend analysis plot
win.graph(height=8,width=5)
par(las=1)
par(mfcol=c(3,1))

plot(Year[26:57],AdjPDI[26:57],type="l",xlab="",ylab="Atlantic PDI, 0x",ylim=c(0,35))
points(Year[26:57],AdjPDI[26:57],pch=21)
lines(Year[26:57],PDI.lm.predict$fit[,1])
lines(Year[26:57],PDI.lm.predict$fit[,2])
lines(Year[26:57],PDI.lm.predict$fit[,3])
abline(v=c(1950,1960,1970,1980,1990,2000),col="grey")
text(1977,28,"(r.sq=0.28, p=0.0019)")

plot(Year[26:57],AdjPDI.flt[26:57],type="l",xlab="",ylab="Atlantic PDI, 1x",ylim=c(0,35))
points(Year[26:57],AdjPDI.flt[26:57],pch=21)
lines(Year[26:57],PDI.flt.lm.predict$fit[,1])
lines(Year[26:57],PDI.flt.lm.predict$fit[,2])
lines(Year[26:57],PDI.flt.lm.predict$fit[,3])
abline(v=c(1950,1960,1970,1980,1990,2000),col="grey")
text(1977,28,"(r.sq=0.43, p=0.0000)")

plot(Year[26:57],AdjPDI.2.flt[26:57],type="l",xlab="",ylab="Atlantic PDI, 2x",ylim=c(0,35))
points(Year[26:57],AdjPDI.2.flt[26:57],pch=21)
lines(Year[26:57],PDI.2.flt.lm.predict$fit[,1])
lines(Year[26:57],PDI.2.flt.lm.predict$fit[,2])
lines(Year[26:57],PDI.2.flt.lm.predict$fit[,3])
abline(v=c(1950,1960,1970,1980,1990,2000),col="grey")
text(1977,28,"(r.sq=0.51, p=0.0000)")