library(fExtremes)
#library(evir)
library(ismev)


########
#gev eloszlasok
########
par(mfrow=c(1,1))
n=seq(-2,6,0.01)
plot(c(n),dgev(n,xi=0.8,mu=0,beta=1), main="Sűrűségfüggvények", ylab="y",xlab="x", type="l", lwd=1.5)

#xi szerepe

lines(c(n),dgev(n,xi=0.6,mu=0,beta=1), col=2, lwd=1.5)
lines(c(n),dgev(n,xi=0.4,mu=0,beta=1), col=3, lwd=1.5)

lines(c(n),dgev(n,xi=0.2,mu=0,beta=1), col=4, lwd=1.5)
lines(c(n),dgev(n,xi=0.2,mu=1,beta=2), col=4, lwd=1.5)

lines(c(n),dgev(n,xi=0,mu=0,beta=1), col=7, lwd=2.5)

legend("topright",c("xi=0.8","xi=0.6","xi=0.4","xi=0.2","xi=0"),lty=c(1,1,1,1,1),col=c(1,2,3,4,7))


lines(c(n),dgev(n,xi=0,mu=0,beta=1), col=7, lwd=2.5)

legend("topright",c("xi=0.8","xi=0.6","xi=0.4","xi=0.2","xi=0"),lty=c(1,1,1,1,1),col=c(1,2,3,4,7))

#nezzuk meg a tobbi parametert is

#szimul
n=c(25,100,400)

j=3
r.gev=rgev(n[j],xi=0.4,mu=0,beta=1) 
becs=gev.fit(r.gev)
gev.diag(becs)

rep=100
ered=array(0,dim=c(rep,3,4))

for (j in 1:3)
{

for (i in 1:rep)
{r.gev=rgev(n[j],xi=0.6,mu=0,beta=1) 
ered[i,j,1:3]=gev.fit(r.gev)$mle
ered[i,j,4]=qgev(.95,mu=ered[i,j,1],beta=ered[i,j,2],xi=ered[i,j,3]) #95%-os VaR

}
}

par(mfrow=c(2,2))
for (i in 1:4)
hist(ered[,3,i])#mik is ezek?

#nezzuk meg mas parameterekre is

#VaR becsles
j=2
par(mfrow=c(1,1))
r.gev=rgev(n[j],xi=0.6,mu=0,beta=1) 
becs=gev.fit(r.gev)
gev.prof(becs,20,xlow=5,xup=24)
 #95%-os. Mi lenne a pontos ertek?

qgev(0.95,xi=0.6,mu=0,beta=1) 
###############
#adatok
################

library(fExtremes)
#library(evir)
library(ismev)



y=read.table("D:\\OKTATAS\\nasdaq_logreturn.txt")
y1=as.matrix(y,nc=1)
yy=-y1 #losses!

# block maxima method
b=blockMaxima(yy,20) #monthly 
hist(b)

#Fitting distribution

#f=gevFit(b,1)

f=gevFit(yy,20) #because that fit is based on the same method

#illustration
xi=f@fit$par.ests[1]
mu=f@fit$par.ests[2]
beta=f@fit$par.ests[3]
h=hist(b,prob=TRUE,ylim=c(0,35))
xfit=seq(min(b),max(b),length=40) 
yfit=dgev(xfit,xi,mu,beta)
lines(xfit,yfit,col="red",ylim=c(min(yfit),max(yfit)))

#goodness of fit
#Diagnostics
zz=gev.fit(b) 
gev.diag(zz)

z=seq(0,1,0.1)

vec=qgev(z,xi=xi,mu=mu,beta=beta)
bb.cut=cut(b,breaks=vec)
table(bb.cut)

f.ex=rep(1/(length(z)-1),times=(length(z)-1))*length(b)

f.obs=table(bb.cut)
X2=sum(((f.obs-f.ex)^2)/f.ex)
df=(length(vec)-1)-3-1
1-pchisq(X2,df)


#VaR calculation
qgev(0.95,xi,mu,beta)
zz=gev.fit(b); #goodness of fit
gev.diag(zz); #Diagnostics


gev.profxi(zz,0,0.61,conf=0.95,nint=100); #profile lik. of shape parameter 
gev.prof(zz,m=60,0.06,0.15,conf=0.95,nint=100); #profile lik. of shape parameter 
qgev(1-1/60,xi,mu,beta); #VaR calculation