#Load necessary packages
library(extRemes)


#Read the data
AuroraMonthlyMaxSeasonal <- read.table("AuroraMonthlyMaxSeasonal.txt")
SeasonalDemand <- read.table("SeasonalDemand.txt")
CovMatGEVAurora <- as.matrix(read.table("CovMatGEVAurora.txt",header=TRUE))
CovMat <- as.matrix(read.table("CovMat.txt"))


#Stationary Seasonal Monthly 
fit.gev6 <- gev.fit( AuroraMonthlyMaxSeasonal[,1])
gev.diag(fit.gev6)

#Stationary GPD
gpd.fitrange( SeasonalDemand[,1], umin=70, umax=115, nint=15)

gpd.fit0 <- gpd.fit( SeasonalDemand[,1], 90,npy=153)
gpd.diag( gpd.fit0)

#Stationary PP
pp.fitrange( SeasonalDemand[,1], umin=70,umax=115,nint=15)

pp.fit0 <- pp.fit( SeasonalDemand[,1], 90,npy=153)
pp.diag( pp.fit0)

#Non-Stationary GEV
fit.gev7 <- gev.fit( AuroraMonthlyMaxSeasonal[,1],ydat=CovMatGEVAurora, sigl=c(1,2,3,4))
gev.diag(fit.gev7)

fit.gev8 <- gev.fit( AuroraMonthlyMaxSeasonal[,1],ydat=CovMatGEVAurora, mul=c(1,2,3,4))
gev.diag(fit.gev8)

fit.gev78 <- gev.fit( xdat= AuroraMonthlyMaxSeasonal[,1],ydat=CovMatGEVAurora, mul=c(1,2,3,4), sigl=c(1,2,3,4))
gev.diag(fit.gev78)

deviancestat( fit.gev6$nllh, fit.gev7$nllh, v=1)
deviancestat( fit.gev6$nllh, fit.gev8$nllh, v=1)
deviancestat( fit.gev6$nllh, fit.gev78$nllh, v=2)
deviancestat( fit.gev7$nllh, fit.gev78$nllh, v=1)
deviancestat( fit.gev8$nllh, fit.gev78$nllh, v=1)

#Non-Stationary GPD
gpd.fit1 <- gpd.fit(SeasonalDemand[,1], 90, npy=153, ydat = CovMat, sigl = c(1,2,3,4), siglink = identity)
gpd.diag( gpd.fit1)

deviancestat( gpd.fit0$nllh, gpd.fit1$nllh, v=1)

#Non-Stationary PP
pp.fit1 <- pp.fit(SeasonalDemand[,1],95,npy=153,ydat=CovMat, sigl = c(1,2,3,4), siglink = identity)
pp.diag( pp.fit1)

pp.fit3 <- pp.fit(SeasonalDemand[,1],95,npy=153,ydat=CovMat, mul=c(1,2,3,4))
pp.diag( pp.fit3)

pp.fit3 <- pp.fit(SeasonalDemand[,1],95,npy=153,ydat=CovMat, mul=c(1,2,3,4))
pp.diag( pp.fit3)

deviancestat( pp.fit0$nllh, pp.fit1$nllh, v=1)
deviancestat( pp.fit0$nllh, pp.fit3$nllh, v=1)
deviancestat( pp.fit0$nllh, pp.fit13$nllh, v=2)
deviancestat( pp.fit1$nllh, pp.fit13$nllh, v=1)
deviancestat( pp.fit3$nllh, pp.fit13$nllh, v=1)