attach(data)
data[1:10,] #prohledneme si prvnich 10 zaznamu

library(fExtremes)
library(evd)

#(1)
A<-data$A   #dale budeme pracovat jen s dennimi srazkami na miste A

#a)
A<-A[1:16071]
plot(A,type="l")

#b)
mrlplot(A,c(1,100)) 

#c)
u<-quantile(A,0.96)[[1]]    
u
plot(A,type="l")
abline(h=u,col=2)

#d)
fpot(A,u)   #odhady parametru GPD rozdeleni
confint(fpot(A,u)) #intervaly spolehlivosti

gpdFit(A, u = u, type = "mle") #jina funkce
gpdFit(A, u = u, type = "pwm") #metoda pravd. vazenych momentu

sigma1=fpot(A,u)$est[[1]]
gamma1=fpot(A,u)$est[[2]]

gamma1  #porovejme s odhadem parametru gamma GEV rozdeleni (0.2322)

N<-length(A)   #puvodni pocet pozorovani
Y<-A[A>u]-u    #excesy
Nu<-length(Y)  #pocet excesu

#e)
hist(Y,prob=T,ylim=c(0,0.11),n=11,main="Histogram excesu")
curve(evd::dgpd(x,0,sigma1,gamma1),col=2,add=T) #s odhadnutou hustotou GPD rozdeleni

plot(density(Y),ylim=c(0,0.11),main="Jadrovy odhad hustoty excesu")
curve(evd::dgpd(x,0,sigma1,gamma1),0,130,col=2,add=T) #s odhadnutou hustotou GPD rozdeleni

library(logKDE)
plot(logdensity(Y),ylim=c(0,0.18),main="Jadrovy odhad hustoty excesu")
curve(evd::dgpd(x,0,sigma1,gamma1),0,130,col=2,add=T) #s odhadnutou hustotou GPD rozdeleni

plot(ecdf(Y),main="Empiricka distribucni funkce excesu")
curve(evd::pgpd(x,0,sigma1,gamma1),col=2,add=TRUE)

# Q-Q plot
i<-1:Nu
q<-evd::qgpd((i-0.5)/Nu,0,sigma1,gamma1)
plot(sort(Y)~q,xlab='Theoretical quantiles',ylab='Empirical quantiles', main='Q-Q plot')
abline(0,1)

#P-P plot
p<-i/(Nu+1)
F<-evd::pgpd(sort(Y),0,sigma1,gamma1)
plot(p~F,xlab='Theoretical cdf',ylab='Empirical cdf', main='P-P plot')
abline(0,1)

#f)
# pravdepodobnost, ze srazkovy uhrn dany den prekroci hodnotu x:
x<-100

p<-(1-evd::pgpd(x-u,0,sigma1,gamma1))*Nu/N
p

#g)
# pravdepodobnost, ze maximum srazkovych uhrnu za dany rok prekroci hodnotu x:
1-(1-p)^365

#h)
# uroven navratu - extremni jev, ktery nastane v prumeru 1x za k dni:
k<-100*365 #pocet dnu odpovidajici priblizne 100 letum

u+evd::qgpd(1-N/(k*Nu),0,sigma1,gamma1)

#i)
# doba navratu - prumerna frekvence extremniho jevu dane urovne:
q<-150
k<-N/((1-evd::pgpd(q-u,0,sigma1,gamma1))*Nu)

k        #frekvence ve dnech
k/365    #frekvence v letech 

plot(fpot(A,u))


#(2)
B=data$B

#a)
B<-B[1:14610]
plot(B,type="l")

#b)
mrlplot(B,c(1,100)) 

#c)
u<-quantile(B,0.96)[[1]]    
u
plot(B,type="l")
abline(h=u,col=2)

#d)
fpot(B,u)   #odhady parametru GPD rozdeleni
confint(fpot(B,u)) #intervaly spolehlivosti

gpdFit(B, u = u, type = "mle") #jina funkce
gpdFit(B, u = u, type = "pwm") #metoda pravd. vazenych momentu

sigma1<-fpot(B,u)$est[[1]]
gamma1<-fpot(B,u)$est[[2]]

N<-length(B)   #puvodni pocet pozorovani
Y<-B[B>u]-u    #excesy
Nu<-length(Y)  #pocet excesu

#e)
hist(Y,prob=T,ylim=c(0,0.08),n=11,main="Histogram excesu")
curve(evd::dgpd(x,0,sigma1,gamma1),col=2,add=T) #s odhadnutou hustotou GPD rozdeleni

plot(density(Y),ylim=c(0,0.08),main="Jadrovy odhad hustoty excesu")
curve(evd::dgpd(x,0,sigma1,gamma1),0,130,col=2,add=T) #s odhadnutou hustotou GPD rozdeleni

library(logKDE)
plot(logdensity(Y),ylim=c(0,0.1),main="Jadrovy odhad hustoty excesu")
curve(evd::dgpd(x,0,sigma1,gamma1),0,130,col=2,add=T) #s odhadnutou hustotou GPD rozdeleni

plot(ecdf(Y),main="Empiricka distribucni funkce excesu")
curve(evd::pgpd(x,0,sigma1,gamma1),col=2,add=TRUE)

# Q-Q plot
i<-1:Nu
q<-evd::qgpd((i-0.5)/Nu,0,sigma1,gamma1)
plot(sort(Y)~q,xlab='Theoretical quantiles',ylab='Empirical quantiles', main='Q-Q plot')
abline(0,1)

#P-P plot
p<-i/(Nu+1)
F<-evd::pgpd(sort(Y),0,sigma1,gamma1)
plot(p~F,xlab='Theoretical cdf',ylab='Empirical cdf', main='P-P plot')
abline(0,1)

#f)
# pravdepodobnost, ze srazkovy uhrn dany den prekroci hodnotu x:
x<-100

p<-(1-evd::pgpd(x-u,0,sigma1,gamma1))*Nu/N
p

#g)
# pravdepodobnost, ze maximum srazkovych uhrnu za dany rok prekroci hodnotu x:
1-(1-p)^365

#h)
# uroven navratu - extremni jev, ktery nastane v prumeru 1x za k dni:
k<-100*365 #pocet dnu odpovidajici priblizne 100 letum

u+sigma1/gamma1*((k*Nu/N)^(gamma1)-1)

#i)
# doba navratu - prumerna frekvence extremniho jevu dane urovne:
q<-150
k<-1/((1-evd::pgpd(q-u,0,sigma1,gamma1))*Nu/N)

k        #frekvence ve dnech
k/365    #frekvence v letech