load(evd,mev,ismev,scales,lubridate,gridExtra,ggplot2,dplyr,tidyr,ggdist,ggpubr,xts)

load("eskrain.RData")

time.seq <- seq(from=min(date(eskrain)), to=max(date(eskrain)), length=149016)
precip_numeric <- as.numeric(eskrain)
esk.rain <- data.frame(date=as.Date(time.seq), precip=precip_numeric)

u<-5
plot_esk <- plot(esk.rain, type="h", ylab="Hourly rainfall (mm)", xlab="Time")
points(esk.rain[esk.rain$precip > u,], col="red", cex=.25, pch=20)
abline(u, 0, col="red")

### start by taking daily maxima
daily.max<-apply(matrix(esk.rain$precip, ncol=24, byrow=T), 1, function(x)max(x))

### then use daily maxima to compute weekly maxima
weekly.max<-apply(matrix(daily.max, ncol=7, byrow=T), 1, function(x)max(x))

### finally, take monthly maxima, with months of 30 days (!)
monthly.max<-apply(matrix(daily.max, ncol=30, byrow=T), 1, function(x)max(x))

fit.weekly<- fgev(weekly.max)
par(mfrow=c(2, 2))
plot(fit.weekly)

par(mfrow=c(1,3))
plot(profile(fit.monthly))

fit.monthly0<- fgev(monthly.max, shape=0)
ratio<-fit.monthly0$dev-fit.monthly$dev
qchisq(0.95,1)
p<-1-pchisq(ratio1,1)