StormMax=read.csv("http://myweb.fsu.edu/jelsner/extspace/extremedatasince1899.csv",T)
A<-unlist(StormMax)
write(A,file="HurricaneData",ncolumns=6,sep=",")
#
StormMaxBasin=subset(StormMax,Region=="Basin")
StormMaxBasin=subset(StormMaxBasin,Yr>1977)

attach(StormMaxBasin)
x<-boxplot(Wmax~as.factor(Yr),plot=F)
boxplot(Wmax~as.factor(Yr),ylim=c(35,175),xlab="Year",ylab="Intensity (kt)")
xx=1:32
abline(lm(x$stats[5,]~xx),col="red")
abline(lm(x$stats[4,]~xx),col="blue")
abline(lm(x$stats[3,]~xx),col="green")
########################################################
con = "http://www.hurricaneclimate.com/storage/chapter-8/LMI.txt"
LMI.df = read.table(con, header = TRUE)
round(head(LMI.df)[c(1, 5:9, 12, 16)], 1)
LMI.df$WmaxS = LMI.df$WmaxS * 0.5144
#
quantile(LMI.df$WmaxS, c(0.25, 0.75))
#
pt = seq(0, 1, 0.01)
qp = quantile(LMI.df$WmaxS, prob = pt)
par(mfrow = c(1, 2), las = 1, mgp = c(2, 0.4, 0), tcl = -0.3)
plot(ecdf(LMI.df$WmaxS), xlab = "Wind speed [m/s]", ylab = "Cumulative distribution", 
    main = "")
abline(v = 40, col = "red")
abline(v = 65, col = "red")
rug(LMI.df$WmaxS)
mtext("a", side = 3, line = 1, adj = 0, cex = 1.1)
plot(pt, qp, type = "l", xlab = "Quantile", ylab = "Wind speed [m/s]", lwd = 2)
mtext("b", side = 3, line = 1, adj = 0, cex = 1.1)
#
boxplot(LMI.df$WmaxS ~ as.factor(LMI.df$SYear))
##########################################################
con = "http://www.hurricaneclimate.com/storage/chapter-8/SST.txt"
SST = read.table(con, header = TRUE)
lg = SST$Year >= 1967
sst.df = data.frame(Yr = SST$Year[lg], sst = SST$Jun[lg])
#
lmisst.df = merge(LMI.df, sst.df, by = "Yr")
head(lmisst.df[c("Yr", "WmaxS", "sst")])
#
brk = quantile(lmisst.df$sst, prob = seq(0, 1, 0.2))
sst.i = cut(lmisst.df$sst, brk, include.lowest = TRUE)
#
con = "http://www.hurricaneclimate.com/storage/chapter-8/SOI.txt"
SOI = read.table(con, header = TRUE)
lg = SOI$Year >= 1967
soi.df = data.frame(Yr = SOI$Year[lg], soi = SOI$Sep[lg])
lmisoi.df = merge(LMI.df, soi.df, by = "Yr")
brk = quantile(lmisoi.df$soi, prob = seq(0, 1, 0.2))
soi.i = cut(lmisoi.df$soi, brk, include.lowest = TRUE)
#
xlabs = c("MB", "B", "N", "A", "MA")
W = lmisst.df$WmaxS
y = boxplot(W ~ sst.i, plot = FALSE)
#
par(mfrow = c(1, 2), las = 1, mgp = c(2, 0.4, 0), tcl = -0.3)
y = boxplot(W ~ sst.i, plot = FALSE)
boxplot(W ~ sst.i, notch = TRUE, names = xlabs, xlab = "June SST pentile", ylab = "Wind speed [m/s]", 
    cex.axis = 0.7)
x = 1:5
abline(lm(y$stats[3, ] ~ x), lwd = 2)
abline(lm(y$stats[4, ] ~ x), col = "red", lwd = 2)
mtext("a", side = 3, line = 1, at = 0, cex = 1.1)
y = boxplot(W ~ soi.i, plot = FALSE)
boxplot(W ~ soi.i, notch = TRUE, names = xlabs, xlab = "September SOI pentile", 
    ylab = "Wind speed [m/s]", cex.axis = 0.7)
abline(lm(y$stats[3, ] ~ x), lwd = 2)
abline(lm(y$stats[4, ] ~ x), col = "red", lwd = 2)
mtext("b", side = 3, line = 1, at = 0, cex = 1.1)
#
require(quantreg)
Year = lmisst.df$Yr
W = lmisst.df$WmaxS
SOI = lmisoi.df$soi
SST = lmisst.df$sst
qrm = rq(W ~ Year + SST + SOI, tau = 0.5)
qrm
#
range(SST)
#
qrm2 = rq(W ~ SST, tau = 0.5)
summary(qrm2)
#
summary(rq(W ~ SST, tau = 0.9), se = "iid")