# preliminary version of analysis of NIDD and Phoenix temperature datsets
# using the extRemes package

# please do not take any of this as definitive!

# NIDD data downloaded from NRFA website (station 27001), see for example

# https://nrfa.ceh.ac.uk/data/search

# https://nrfa.ceh.ac.uk/data/station/info/27001



## Nidd analysis - annual maxima
X=read.delim('C:/Users/rsmith/jan16/UNC/STOR834/Rcode/27001.am.txt')
day=as.numeric(substr(X[,1],1,2))
mon=substr(X[,1],4,6)
year=as.numeric(substr(X[,1],8,11))
flow=as.numeric(substr(X[,1],13,21))
mon[mon=='Jan']=1
mon[mon=='Feb']=2
mon[mon=='Mar']=3
mon[mon=='Apr']=4
mon[mon=='May']=5
mon[mon=='Jun']=6
mon[mon=='Jul']=7
mon[mon=='Aug']=8
mon[mon=='Sep']=9
mon[mon=='Oct']=10
mon[mon=='Nov']=11
mon[mon=='Dec']=12
mon=as.numeric(mon)
yr=(1:84)
NIDDAM=data.frame(year,mon,day,yr,flow)

## Nidd analysis - peaks over threshold
X=read.delim('C:/Users/rsmith/jan16/UNC/STOR834/Rcode/27001.pt1.txt')
day=as.numeric(substr(X[,1],1,2))
mon=substr(X[,1],4,6)
year=as.numeric(substr(X[,1],8,11))
flow=as.numeric(substr(X[,1],13,21))

day1=(year-1934)*365
for(jj in 1:21){day1[year>1932+jj*4]=day1[year>1932+jj*4]+1}
day1[mon=='Feb']=day1[mon=='Feb']+31
day1[mon=='Mar']=day1[mon=='Mar']+59
day1[mon=='Apr']=day1[mon=='Apr']+90
day1[mon=='May']=day1[mon=='May']+120
day1[mon=='Jun']=day1[mon=='Jun']+151
day1[mon=='Jul']=day1[mon=='Jul']+181
day1[mon=='Aug']=day1[mon=='Aug']+212
day1[mon=='Sep']=day1[mon=='Sep']+243
day1[mon=='Oct']=day1[mon=='Oct']+273
day1[mon=='Nov']=day1[mon=='Nov']+304
day1[mon=='Dec']=day1[mon=='Dec']+334
leap=4*floor(year/4)==year
day1[mon!='Jan'&mon!='Feb'&leap]=day1[mon!='Jan'&mon!='Feb'&leap]+1
day1=day1+day
day1=day1-312
NIDDPT=data.frame(day1,flow)


# Initial runs with extRemes package

library(extRemes)
f1=fevd(flow,NIDDAM,type='GEV')
ci.fevd(f1,return.period=500,method='boot')
ci.fevd(f1,return.period=500,method='proflik')
ci.fevd(f1,return.period=500,method='normal')


f2=fevd(flow,NIDDAM,location.fun=~yr,type='GEV')

f3=fevd(flow,NIDDAM,location.fun=~yr+I(yr^2),type='GEV')

f4=fevd(flow,NIDDAM,location.fun=~yr,scale.fun=~yr,type='GEV')


ci(f2,return.period=500,method='normal')

# POT analysis of NISS data

mrlplot(NIDDPT$flow)

f5=fevd(flow,NIDDPT,threshold=70,type='GP')
f6=fevd(flow,NIDDPT,threshold=100,type='GP')

ci.fevd(f5,return.period=100,method='boot')
ci.fevd(f5,return.period=100,method='proflik')
ci.fevd(f5,return.period=100,method='normal')

ci.fevd(f6,return.period=100,method='boot')
ci.fevd(f6,return.period=100,method='proflik')
ci.fevd(f6,return.period=100,method='normal')





# Use Phoenix temperature data from package
data(HEAT)
f1=fevd(Tmax,HEAT,type='GEV')
ci.fevd(f1,return.period=500,method='boot')
ci.fevd(f1,return.period=500,method='proflik')
ci.fevd(f1,return.period=500,method='normal')


f2=fevd(Tmax,HEAT,location.fun=~Year,type='GEV')


f3=fevd(Tmax,HEAT,location.fun=~Year,scale.fun=~Year,type='GEV')


HEAT$negtmin=-HEAT$Tmin



f4=fevd(negtmin,HEAT,type='GEV')