#########################################################################
###          Extremwerttheorie  -  Übungsblatt 5	                ###
###          SS 2017		  Dr. Jürgen Kampf       	    	    ###
#########################################################################


### Aufgabe 2
#############

# Teil a)
# Laden der Pakets und Verfügbarmachen des Datensatzes
library(evd)
data(sealevel)

# Übersichtsplots (nicht verlangt)
plot(1912:1992, sealevel$dover, xlab="Jahr", ylab="Maximale Höhe", pch=16)
plot(ecdf(sealevel$dover))

# Anwenden des Schätzers
fgev(sealevel$dover)
# Unter estimates sehen wir die Resultate des Maximum-Likelihood-Schätzers
# gamma = shape = -0.02107
# mu = loc = 3.59252 
# sigma = scale = 0.20195

# Teil b)
# Entfernen der fehlenden Werte
sea_dover=sealevel$dover[!is.na(sealevel$dover)]

# Log-Likelihood-Funktion Variante 1: Von Hand

LL=function(gamma,mu,sigma){
	if(sigma<0){return(-1000)}	# Unzulässiges Argument => gebe sehr niedrigen Wert zurück
	llsum = 0;				# Summationsvariable für den Wert der Log-Likelihhod-Funktion
	for(z in sea_dover){
	if(gamma==0){
		# Gumbel
		llsum= llsum - log(sigma) - (z-mu)/sigma  - exp(-(z-mu)/sigma)
		}
	else{
		# Frechet/Weibull
		if( 1+(gamma*(z-mu)/sigma) < 0 ){return(-1000)} 
						# Unzulässiges Argument => gebe sehr niedrigen Wert zurück

		else{ 
		llsum  =  llsum  - log(sigma) - (1/gamma+1)*log(1+gamma*(z-mu)/sigma) - (1+gamma*(z-mu)/sigma)^(-1/gamma) 
		}}
	}# for z
	return(llsum)
	}


# Log-Likelihood-Funtion Variante 2: Mit Hilfe des evd-Pakets

LL=function(gamma,mu,sigma){
	if(sigma<0){return(-1000)}	# Unzulässiges Argument => gebe sehr niedrigen Wert zurück
	llsum = 0;				# Summationsvariable für den Wert der Log-Likelihhod-Funktion
	for(z in sea_dover){
		# Frechet/Weibull
		if( 1+(gamma*(z-mu)/sigma) < 0 ){return(-1000)} 
						# Unzulässiges Argument => gebe sehr niedrigen Wert zurück

		else{ 
		llsum  =  llsum  + log(dgev(z, shape=gamma, loc=mu, scale= sigma)) 
		}
	}# for z
	return(llsum)
	}


# Vorbereitende Schritte zur Optimierung
# Die Funktion darf nur ein Argument haben; außerdem muss ihr Negatives betrachtet werden, da
# optim() das Minimum bestimmt, wir aber das Maximum suchen

NLL = function(args){
	if( length(args)!=3){
		print("In NLL: Falsche Länge des Argumentvektors");
		return(NA);
		}
	return( -LL( gamma=args[1], mu=args[2], sigma=args[3]) );
	}

# Optimierung
optim(par=c(0, 2, 1), fn=NLL)
MLE=optim(par=c(0, 2, 1), fn=NLL)$par


# Teil c)
# Plot der Log-Likelihood-Funtion in Abhängigkeit von mu
mu_seq=(201:400)/100
ll_seq=rep(0,200)
for(i in (1:200)){
	ll_seq[i]=LL(MLE[1], mu_seq[i], MLE[3])
	}
plot(mu_seq,ll_seq,type="l",xlab="mu",ylab="l")

# Plot der Log-Likelihood-Funktion in Abhängigkeit von sigma 
sigma_seq=(1:400)/1000
ll_seq=rep(0,400)
for(i in (1:400)){
	ll_seq[i]=LL(MLE[1], MLE[2], sigma_seq[i])
	}
plot(sigma_seq,ll_seq,type="l", ylim=c(-200,0), xlab="sigma", ylab="l")


# Teil d)
# Der QQ-Plot
len=length(sea_dover)
p=(1:len)/(len+1)
plot(qgev(p,loc=MLE[2],scale=MLE[3],shape=MLE[1]),sort(sea_dover), xlab="Theoretische Quantile",
		ylab="Emprische Quantile")
abline(a=0,b=1)
# Die Punkte liegen ziemlich exakt auf der 1. Winkelhalbierenden.
# => Die Wasserstände lassen sich gut durch eine GEV-Verteilung beschreiben.