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

#(1)
A=data$A   #zacneme s dennimi srazkami na miste A

#a)
plot(A,type="l")

par(mfrow=c(2,2))
plot(A[1:365], main="1961",ylim=c(0,60),type="l")
plot(A[366:730], main="1962",ylim=c(0,60),type="l")
plot(A[731:1095], main="1963",ylim=c(0,60),type="l")
plot(A[1096:1461], main="1964",ylim=c(0,60),type="l")

hist(A,n=50,xlab="denni srazky",freq=FALSE)

#b)
n=length(Rok)
a=min(Rok)
b=max(Rok)
X=rep(-1,b-a+1)

for (i in 1:n){
r=Rok[i]
if (A[i]>X[r-a+1]) X[r-a+1]=A[i]}
X=X[-45] #odstranime posledni pozorovani odpovidajici r. 2005 (jen par mesicu)
X

par(mfrow=c(1,2))
rok=(1961:2004)
N=length(rok)
plot(X~rok,main="Rocni maximalni denni srazky")
plot(X~rok,type="l",main="Rocni maximalni denni srazky")

#c)
library(evd)
fgev(X)  #odhad parametru GEV rozdeleni metodou max. verohodnosti
confint(fgev(X))  #intervaly spolehlivosti pro parametry GEV rozdeleni

gamma=fgev(X)$est[[3]]
mu=fgev(X)$est[[1]]
sigma=fgev(X)$est[[2]]

# srovnani s jinou funkci + metoda PWM
library(fExtremes)
gevFit(X, type = "mle")
gevFit(X, type = "pwm")


#d)
hist(X,prob=T,ylim=c(0,0.04),n=11,main="Histogram blokovych maxim") #histogram blokovych maxim
curve(evd::dgev(x,mu,sigma,gamma),0,130,col=2,add=T) #s odhadnutou hustotou GEV rozdeleni

plot(density(X),ylim=c(0,0.04),main="Jadrovy odhad hustoty blokovych maxim")
curve(evd::dgev(x,mu,sigma,gamma),0,130,col=2,add=T) #s odhadnutou hustotou GEV rozdeleni

plot(ecdf(X),main="Empiricka distribucni funkce blokovych maxim")
curve(evd::pgev(x,mu,sigma,gamma),col=2,add=TRUE)

# Q-Q plot
i=1:N
q=evd::qgev((i-0.5)/N,mu,sigma,gamma)
plot(sort(X)~q,xlab='Theoretical quantiles',ylab='Empirical quantiles', main='Q-Q plot')
abline(0,1)

#P-P plot
p=i/(N+1)
F=evd::pgev(sort(X),mu,sigma,gamma)
plot(p~F,xlab='Theoretical cdf',ylab='Empirical cdf', main='P-P plot')
abline(0,1)

#e)
# pravdepodobnost, ze maximum srazkovych uhrnu za dany rok prekroci hodnotu x:
x=100
q=1-evd::pgev(x,mu,sigma,gamma)
q

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

1-evd::pgev(x,mu,sigma,gamma)^(1/365) # presny vzorec
(1+gamma*(x-mu)/sigma)^(-1/gamma)/365 # aproximace

#g)
# uroven navratu - extremni jev, ktery nastane v prumeru 1x za k let:
k=100 #stolete srazky
evd::qgev(1-1/k,mu,sigma,gamma)

#h)
# doba navratu - prumerna frekvence extremniho jevu dane urovne (v letech):
q=150 #denni srazky vyssi nez 150mm (v letech)
1/(1-evd::pgev(q,mu,sigma,gamma))

plot(fgev(X))


####################################################################
####################################################################
#(2)
B=data$B   #a nyni denni srazky na miste B

#a)
plot(B,type="l")

par(mfrow=c(2,2))
plot(B[1:365], main="1961",ylim=c(0,80),type="l")
plot(B[366:730], main="1962",ylim=c(0,80),type="l")
plot(B[731:1095], main="1963",ylim=c(0,80),type="l")
plot(B[1096:1461], main="1964",ylim=c(0,80),type="l")

hist(B,n=50,xlab="denni srazky",freq=FALSE)

#b)
n=length(Rok)
a=min(Rok)
b=max(Rok)
X=rep(-1,b-a+1)

for (i in 1:n){
r=Rok[i]
if (B[i]>X[r-a+1]) X[r-a+1]=B[i]}
X=X[-c(41:45)] #odstranime poslednich 5 pozorovani odpovidajici r. 2001 - 2005
X

par(mfrow=c(1,2))
rok=(1961:2000)
N=length(rok)
plot(X~rok,main="Rocni maximalni denni srazky")
plot(X~rok,type="l",main="Rocni maximalni denni srazky")

#c)
fgev(X)  #odhad parametru GEV rozdeleni metodou max. verohodnosti
confint(fgev(X))  #intervaly spolehlivosti pro parametry GEV rozdeleni

gamma=fgev(X)$est[[3]]
mu=fgev(X)$est[[1]]
sigma=fgev(X)$est[[2]]

#d)
hist(X,prob=T,ylim=c(0,0.025),n=11,main="Histogram blokovych maxim") #histogram blokovych maxim
curve(evd::dgev(x,mu,sigma,gamma),0,130,col=2,add=T) #s odhadnutou hustotou GEV rozdeleni

plot(density(X),ylim=c(0,0.025),main="Jadrovy odhad hustoty blokovych maxim")
curve(evd::dgev(x,mu,sigma,gamma),0,130,col=2,add=T) #s odhadnutou hustotou GEV rozdeleni

plot(ecdf(X),main="Empiricka distribucni funkce blokovych maxim")
curve(evd::pgev(x,mu,sigma,gamma),col=2,add=TRUE)

# Q-Q plot
i=1:N
q=evd::qgev((i-0.5)/N,mu,sigma,gamma)
plot(sort(X)~q,xlab='Theoretical quantiles',ylab='Empirical quantiles', main='Q-Q plot')
abline(0,1)

#P-P plot
p=i/(N+1)
F=evd::pgev(sort(X),mu,sigma,gamma)
plot(p~F,xlab='Theoretical cdf',ylab='Empirical cdf', main='P-P plot')
abline(0,1)

#e)
# pravdepodobnost, ze maximum srazkovych uhrnu za dany rok prekroci hodnotu x:
x=100
q=1-evd::pgev(x,mu,sigma,gamma)
q

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

1-evd::pgev(x,mu,sigma,gamma)^(1/365) # presny vzorec
(1+gamma*(x-mu)/sigma)^(-1/gamma)/365 # aproximace

#g)
# uroven navratu - extremni jev, ktery nastane v prumeru 1x za k let:
k=100 #stolete srazky
evd::qgev(1-1/k,mu,sigma,gamma)

#h)
# doba navratu - prumerna frekvence extremniho jevu dane urovne:
q=150 #denni srazky vyssi nez 150mm
1/(1-evd::pgev(q,mu,sigma,gamma))

plot(fgev(X))