# Chicago Marathon analysis

# First load rlarg.test function from Rvenice-rev1 file ("Revised R code to include tests of fit")

# load data

D=read.csv('https://rls.sites.oasis.unc.edu/s834-2023/Data/ChicagoMarathonData.csv',header=T)

# subtract marathon times from 200 to focus on maxima
V=200-D[1:25,2:26]

library(ismev)

# rr is the chosen value of r - coud be 5, 10, 15, 20, 25
rr=25

mod0=rlarg.fit(V, r = rr, ydat = NULL, mul = NULL, sigl = NULL,
  shl = NULL, mulink = identity, siglink = identity, shlink = identity,
  muinit = mean(V[,1]), siginit = sd(V[,1]), shinit = 0.001, show = TRUE,
  method = "Nelder-Mead", maxit = 10000)

mod1=rlarg.fit(V, r = rr, ydat = matrix(1:25,ncol=1), mul = 1, sigl = NULL,
  shl = NULL, mulink = identity, siglink = identity, shlink = identity,
  muinit = c(mean(V[,1]),0), siginit = sd(V[,1]), shinit = 0.01, show = TRUE,
  method = "Nelder-Mead", maxit = 10000)

mod2=rlarg.fit(V, r = rr, ydat = matrix(c(1:25,rep(0,19),1:6),ncol=2), mul = c(1,2), sigl = NULL,
  shl = NULL, mulink = identity, siglink = identity, shlink = identity,
  muinit = c(mod1$mle[1],mod1$mle[2],0), siginit = mod1$mle[3], shinit = mod1$mle[4], show = TRUE,
  method = "Nelder-Mead", maxit = 10000)

mod3=rlarg.fit(V, r = rr, ydat = matrix(c(1:25,rep(0,18),1:7),ncol=2), mul = c(1,2), sigl = NULL,
  shl = NULL, mulink = identity, siglink = identity, shlink = identity,
  muinit = c(mod1$mle[1],mod1$mle[2],0), siginit = mod1$mle[3], shinit = mod1$mle[4], show = TRUE,
  method = "Nelder-Mead", maxit = 10000)

mod4=rlarg.fit(V, r = rr, ydat = matrix(c(1:25,rep(0,17),1:8),ncol=2), mul = c(1,2), sigl = NULL,
  shl = NULL, mulink = identity, siglink = identity, shlink = identity,
  muinit = c(mod1$mle[1],mod1$mle[2],0), siginit = mod1$mle[3], shinit = mod1$mle[4], show = TRUE,
  method = "Nelder-Mead", maxit = 10000)

mod5=rlarg.fit(V, r = rr, ydat = matrix(c(1:25,rep(0,16),1:9),ncol=2), mul = c(1,2), sigl = NULL,
  shl = NULL, mulink = identity, siglink = identity, shlink = identity,
  muinit = c(mod1$mle[1],mod1$mle[2],0), siginit = mod1$mle[3], shinit = mod1$mle[4], show = TRUE,
  method = "Nelder-Mead", maxit = 10000)

mod6=rlarg.fit(V, r = rr, ydat = matrix(c(1:25,rep(0,15),1:10),ncol=2), mul = c(1,2), sigl = NULL,
  shl = NULL, mulink = identity, siglink = identity, shlink = identity,
  muinit = c(mod1$mle[1],mod1$mle[2],0), siginit = mod1$mle[3], shinit = mod1$mle[4], show = TRUE,
  method = "Nelder-Mead", maxit = 10000)

mod7=rlarg.fit(V, r = rr, ydat = matrix(c(1:25,rep(0,14),1:11),ncol=2), mul = c(1,2), sigl = NULL,
  shl = NULL, mulink = identity, siglink = identity, shlink = identity,
  muinit = c(mod1$mle[1],mod1$mle[2],0), siginit = mod1$mle[3], shinit = mod1$mle[4], show = TRUE,
  method = "Nelder-Mead", maxit = 10000)

mod8=rlarg.fit(V, r = rr, ydat = matrix(c(1:25,rep(0,13),1:12),ncol=2), mul = c(1,2), sigl = NULL,
  shl = NULL, mulink = identity, siglink = identity, shlink = identity,
  muinit = c(mod1$mle[1],mod1$mle[2],0), siginit = mod1$mle[3], shinit = mod1$mle[4], show = TRUE,
  method = "Nelder-Mead", maxit = 10000)

mod9=rlarg.fit(V, r = rr, ydat = matrix(c(1:25,rep(0,12),1:13),ncol=2), mul = c(1,2), sigl = NULL,
  shl = NULL, mulink = identity, siglink = identity, shlink = identity,
  muinit = c(mod1$mle[1],mod1$mle[2],0), siginit = mod1$mle[3], shinit = mod1$mle[4], show = TRUE,
  method = "Nelder-Mead", maxit = 10000)

mod10=rlarg.fit(V, r = rr, ydat = matrix(c(1:25,rep(0,11),1:14),ncol=2), mul = c(1,2), sigl = NULL,
  shl = NULL, mulink = identity, siglink = identity, shlink = identity,
  muinit = c(mod1$mle[1],mod1$mle[2],0), siginit = mod1$mle[3], shinit = mod1$mle[4], show = TRUE,
  method = "Nelder-Mead", maxit = 10000)


print(c(mod0$nllh,mod1$nllh,mod2$nllh,mod3$nllh,mod4$nllh,mod5$nllh,mod6$nllh,mod7$nllh,mod8$nllh,mod9$nllh,mod10$nllh))
# [1] 394.1971 384.5971 376.5111 375.7590 375.3602 375.6681 375.1977 374.7875
# [9] 374.6898 374.8736 376.0200
# based on r=25:
# [1] 692.8916 678.3027 669.7693 669.3257 668.9826 668.9884 668.6297 668.3710
# [9] 668.3018 668.0468 668.4783

# test fit of most important models

T1=Sys.time()
gof1=rlarg.test(V, r = 25, ydat = matrix(1:25,ncol=1), mul = 1, sigl = NULL,
  shl = NULL, mulink = identity, siglink = identity, shlink = identity,
  muinit = c(mean(V[,1]),0), siginit = sd(V[,1]), shinit = 0.01, show = FALSE,
  method = "Nelder-Mead", maxit = 10000, nsim=1000)
T2=Sys.time()
print(T2-T1)
gof1

T1=Sys.time()
gof8=rlarg.test(V, r = 25, ydat = matrix(c(1:25,rep(0,13),1:12),ncol=2), mul = c(1,2), sigl = NULL,
  shl = NULL, mulink = identity, siglink = identity, shlink = identity,
  muinit = c(mod1$mle[1],mod1$mle[2],0), siginit = mod1$mle[3], shinit = mod1$mle[4], show = FALSE,
  method = "Nelder-Mead", maxit = 10000, nsim=1000)
T2=Sys.time()
print(T2-T1)
gof8

T1=Sys.time()
gof2=rlarg.test(V, r = 25, ydat = matrix(c(1:25,rep(0,19),1:6),ncol=2), mul = c(1,2), sigl = NULL,
  shl = NULL, mulink = identity, siglink = identity, shlink = identity,
  muinit = c(mod1$mle[1],mod1$mle[2],0), siginit = mod1$mle[3], shinit = mod1$mle[4], show = FALSE,
  method = "Nelder-Mead", maxit = 10000, nsim=1000)
T2=Sys.time()
print(T2-T1)
gof2