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

A=data$A   #dale budeme pracovat jen s dennimi srazkami na miste A
plot(A)

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

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

#pomocna procedura pro vytvoreni blokovych maxim (delka bloku je 1 rok)
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)
plot(X~rok,main="Rocni maximalni denni srazky")
plot(X~rok,type="l",main="Rocni maximalni denni srazky")


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]]

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

# Q-Q plot
i=seq(1:44)
h=i/(44+1)
TQ=qgev(h,mu,sigma,gamma) #teoreticke kvantily
EQ=sort(X)                #empiricke kvantily 

plot(EQ~TQ,xlim=c(20,125),ylim=c(20,125),main="Q-Q-plot")
abline(0,1)

# P-P plot
EDF=(1:44)/(44+1)     #empiricka distribucni funkce
Q=sort(X)
TDF=pgev(Q,mu,sigma,gamma) #teoreticka distribucni funkce

plot(EDF~TDF,xlim=c(0,1),ylim=c(0,1),main="P-P-plot")
abline(0,1)

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


# pravdepodobnost, ze srazkovy uhrn dany den prekroci hodnotu x:
x=100
p=(1+gamma*(x-mu)/sigma)^(-1/gamma)/365 # aproximace
p
1-(1-q)^(1/365) # presny vzorec


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

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


plot(fgev(X))

################################################
metoda POT
#################################################
mrlplot(A,c(1,100)) #graf pro urceni vhodneho prahu
quantile(A,0.96)    #zkusime radsi 96% kvantil
threshold=12        #zvolme hodnotu prahu u=12
fpot(A,threshold)   #odhady parametru zobecneneho paretova rozdeleni

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

gamma   #porovnani odhadu parametru polohy GEV rozdeleni
gamma1  #porovnani odhadu parametru polohy zob. Paretova rozdeleni


#porovnejme empiricke o odhadnute rozdeleni excesu
hist((A-threshold)[A>threshold],prob=T,ylim=c(0,0.11))
curve(dgpd(x,0,sigma1,gamma1),0,130,col=2,add=T)

# pravdepodobnost, ze srazkovy uhrn dany den prekroci hodnotu x:
x=100
pomer=sum(A>threshold)/length(A)
p=(1-pgpd(x-threshold,0,sigma1,gamma1))*pomer
p

# pravdepodobnost, ze maximum srazkovych uhrnu za dany rok prekroci hodnotu x:
x=100
p=(1-pgpd(x-threshold,0,sigma1,gamma1))*sum(A>threshold)/length(A)

1-(1-p)^365

# uroven navratu - extremni jev, ktery nastane v prumeru 1x za k dni:
k=100*365 #priblizne 100 let
alpha=1/k
threshold+sigma1/gamma1*((alpha/pomer)^(-gamma1)-1)

# doba navratu - prumerna frekvence extremniho jevu dane urovne:
uroven=100
p=(1-pgpd(uroven-threshold,0,sigma1,gamma1))*pomer
1/p        #frekvence ve dnech
1/(365*p)  #frekvence v letech