attach(data)
data[1:10,]

A=data$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")

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

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


library(evd)
fgev(X)
confint(fgev(X))

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)
curve(dgev(x,mu,sigma,gamma),0,130,col=2,add=T)

# q-q plot
i=seq(1:44)
h=i/(44+1)

TQ=qgev(h,mu,sigma,gamma)
EQ=sort(X)
plot(EQ~TQ,xlim=c(20,125),ylim=c(20,125),main="Q-Q-plot")
abline(0,1)

# p-p plot
i=rep(1,44)
EDF=cumsum(i)/(44+1)

Q=sort(X)
TDF=pgev(Q,mu,sigma,gamma)
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
qgev(1-1/k,mu,sigma,gamma)

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


plot(fgev(X))

################################################
metoda POT
#################################################
mrlplot(A,c(1,100))
quantile(A,0.96)
threshold=12
fpot(A,threshold)

sigma=fpot(A,threshold)$est[1]
gamma=fpot(A,threshold)$est[2]

hist((A-threshold)[A>threshold],prob=T,ylim=c(0,0.11))
curve(dgpd(x,0,sigma,gamma),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,sigma,gamma))*pomer
p

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

1-(1-p)^365

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

# doba navratu - prumerna frekvence extremniho jevu dane urovne:
u=100
p=(1-pgpd(u-threshold,0,sigma,gamma))*pomer
1/p/365