#####################################################################
# SS 2015, Stochastische Rekursionsgleichungen, Blatt 6, Aufg 5 #
#####################################################################

# Funktion, die den Pfad bis zum Zeitpunkt n eines AR(1) Prozesses mit ARCH(1) Errors simuliert. Dabei:
# - Innovationsprozess ist iid und Standardnormalverteilt
# - Startwert zum Zeitpunkt 0 ist 0
sim<-function(n, alpha, beta, lambda) {
	X<-numeric(n) # Vektor der Laenge n, der nur aus 0'en besteht
	eps<-rnorm(n) # Innovationsprozess
	X[1]<-0 # Startwert von X
	for(i in 2:n) X[i] <- alpha*X[i-1] + eps[i] * sqrt(beta + lambda * X[i-1]^2) # rekursive Berechnung von X_n
	return(X)
}

# Funktion, die die Regression des Datensatzes sim.data durchfuehrt.
# - min.t ist das (moeglichst grosse) t, ab dem die Stuetzstellen gesetzt werden
# - anz.stzstellen ist die Anzahl dieser Stuetzstellen
# Die Stellen werden dabei so gewaehlt, dass das Intervall [min.t, t*] in gleichgrosse Teile zerlegt wird,
# dabei ist t* ein Wert nahe (bis auf 0.001) dem letzten Sprung der empirischen Verteilungsfunktion F der Daten
# (das ist noetig, damit log(1-F(t*)) endlich bleibt).
sim.regression<-function(sim.data, min.t, anz.stzstellen) {
	n<-length(sim.data)
	emp.Vert<-ecdf(sim.data)
	emp.Vert.spruenge<-knots(emp.Vert) # das liefert die Sprungstellen der emp. Verteilungsfunktion
	stzstellen<-seq(min.t,emp.Vert.spruenge[n]-0.001,length.out=anz.stzstellen) # Vektor mit den Stuetzstellen
	plot(log(stzstellen),log(1-emp.Vert(stzstellen)),xlab="log(Zeit)",ylab="log(emp. Tails)")
	model<-lm(log(1-emp.Vert(stzstellen)) ~ log(stzstellen)) # an lineares Modell fitten
	abline(model,col=2)	# Regressionsgrade
	kappa<--model$coefficients[2]; names(kappa)<-c("kappa")
	C<-exp(model$coefficients[1]); names(C)<-c("C")
	legend("topright",legend=paste(c("kappa =","C ="), round(c(kappa,C),3)), box.lty=0,adj=1)
	return(c(kappa,C))
}

# Parameter, mit denen die Funktion aufgerufen werden soll:
n<-10000 # Anzahl Schritte
beta<-0.9
alpha<-0.4
lambda<-0.8 
# Damit wird das (theoretische) kappa = 2.1 [siehe Tabelle 3 in http://www.jstor.org/stable/2699915]

# Um das Ergebnis der Ausgabe reproduzieren zu koennen setzen wir einen 
#  (beliebigen) Startwert fuer die Erzeugung der (Pseudo-)Zufallszahlen mit Hilfe von set.seed()
set.seed(1)

data<-sim(n=n,alpha=alpha,beta=beta,lambda=lambda) # Datensatz erzeugen
sim.regression(data,5,100) # Regression. Liefert hier kappa = 2.3 und C = 0.9