#################################################
### Dateimanagement
#################################################

## ... selbst überlegen, wie es übersichtlich und effizient geht,
## aber als Tipp ...

# Oft verwendete Funktionen (Schätzer, Daten einlesen, ...)
# nicht zu Beginn der Hauptcodes (der Datenanalyse...)
# stellen und in mehreren Files wiederholen, sondern auslagern
# und mit source() einlesen

# Funktionscodes mit Hinweisen versehen (Datum, Funktioniert Code?,
# evtl. Probleme, Verwendung)

# Alternativ eigenes Package erstellen (Tutorial)


#################################################
### Übersichtlicher Code
#################################################
install.packages("AER",repos="http://cran.r-project.org")
install.packages("AER")
library(AER)

### Was macht der folgende Code?

a05 <- 0.05
a1 <- 1
a2 <- 1
a3 <- 1
a4 <- 1
a100 <- 100
a1000 <- 1000

f <- function(a05,a1,a2,a3,a4,a100,a1000) { fi <- function(b05,b1,b2,b3,b4,b100,b123){
d <- list();d$x <- rnorm(b100); set.seed(b123);d$y=d$x*b1+rt(b100,df=b4)*exp(b2+b3*d$x);
 d=as.data.frame(d); library(AER);h <- paste("d$x=",b1,sep="");(linearHypothesis(lm(d$y~d$x),h)[[6]][2]<b05)}
a <- numeric();for (i in 1:a1000) a[i] <- fi(a05,a1,a2,a3,a4,a100,i);mean(a)}
f(0.05,1,1,1,4,100,1000)

### Aufgabe J1: 
# Versuchen Sie, den Code zu verstehen. Bringen Sie ihn dabei in eine Form,
# in der dies einfacher ist. (Code entzerren, Strukturieren, Variablennamen ändern,
# Kommentieren)


#################################################
### Zeitoptimierter Code
#################################################

# Beachte: Erst überlegen, ob sich Optimierung lohnt!
# -> auch Programmierzeit berücksichtigen usw.

### Profiling: Welche aufgerufenen Funktionen bilden den "Flaschenhals"?

Rprof(tmp <- tempfile(), interval = 0.1)
f(0.05,1,1,1,4,100,1000)
Rprof(NULL)
summaryRprof(tmp)


### Prüfen der Laufzeit mit system.time()

system.time(f(0.05,1,1,1,4,100,1000))

### Verringerung der Laufzeit:

### Merke 1: Objekte mit richtiger Größe initialisieren! Kostenaufwand immens!!!!
### Merke 2: Berechnungen aus der Schleife ziehen! Alles unnötige aus der Schleife ziehen!
### Merke 3: Schleifen vermeiden (WO EINFACH MÖGLICH UND SINNVOLL!!!) Ausweg: APPLY
### Merke 4: Unnötige Berechnungen ganz vermeiden
### Merke 5: Implementierte Funktionen verwenden (da schon zeitoptimiert da in C)

### Aufgabe J2: 
# Optimieren Sie die Funktion g schrittweise anhand der Prinzipien 1-5

g <- function(n){
  a <- numeric()
  sa <- 0
  for (i in 1:n){
   a[i] <- sin(i*2*pi)*100
   sa <- sa + a[i]/n
   }
return(sa)
}

### WICHTIG:  
### Merke:
# Erzeugung AR(1) Prozess 


n <- 100000
e <- rnorm(n)
alpha <- 0.9

system.time({x <- e 
             for (i in 2:length(e)) x[i] <- alpha*x[i-1]+e[i]})

system.time(filter(e,alpha,"recursive"))
# Fazit: IMMER MIT FILTER!!!

###
# Regressionsschätzung mit vielen Regressoren

n <- 1500
k <- 1000

X <- matrix(rnorm(n*k),n,k)
y <- rnorm(n)

system.time(beta_hat <- solve(t(X)%*%X)%*%t(X)%*%y)

system.time(beta_hat <- lm(y~-1+X)$coef)

system.time(XtX<-solve(t(X)%*%X))
# Fazit: IMMER MIT MODELL!!!


### Oftmal übersichtlicher und/oder schneller als Schleife:
### apply(), ..., Vectorize()
??"rejections.ttest"
# Schleife
alpha1.seq <- seq(0,1,0.1)
(rej.freq <- seq(along.with=alpha1.seq))

for (i in rej.freq) rej.freq[i] <- t.test(alpha_1=alpha1.seq[i])[1]
rej.freq


# Vectorize
rejections.ttest.V <- Vectorize(rejections.ttest,vectorize.args="alpha_1")

r <- rejections.ttest.V(alpha_1 = seq(0,1,0.1),
  df_error = 1000, n = 250, Reps = 1000)
r

plot(seq(0,1,0.1),r,
  xlab="alpha_1",ylab="Rejection Frequency",
  main="Rejection Probability of True H0 with Heteroskedasticity")
abline(h=5)


### Besonders, wenn Code auch später noch verstanden werden soll
### ODER auch andere ihn verstehen sollen!

# Braucht AER Paket für linear.hypothesis()
library(AER)

# Funktion, die t-Test in linearem Regressionsmodell simuliert
# H_0: wahrer Wert von beta
# t-Verteilter Fehler UND Heteroskedastie (falls alpha_1!=0)
#
# Argumente:
#   sig.level         - Signifikanzniveau
#   beta              - Steigungsparameter
#   alpha_0, alpha_1  - Parameter der exponentiellen Varianzfunktion
#   df_error          - Freiheitsgrade der t-Verteilung bei Simulation von u
#   n                 - Stichprobengröße
#   seed              - Zufallszahlen-Seed
#
# Funktionswert:
#   TRUE für ablehnen
#   FALSE für nicht ablehnen


sim.ttest <- function(sig.level,beta,alpha_0,alpha_1,df_error,n,seed)
{
  set.seed(seed)
  # Generiere Regressor
  x <- rnorm(n)
  # Generiere Fehler: t-Verteilte ZV * Std.Abw (Exponentielle Funktion von x)
  u <- rt(n,df=df_error) * exp(alpha_0+alpha_1*x)
  # Erzeuge Regressanden
  y <- x*beta + u
  
  # Formuliere Nullhypothese
  H0 <- paste("x=",beta,sep="")
  # Schätze Modell mit OLS
  estimated.reg <- lm(y~x)
  # Führe t-Test durch
  ttest <- linear.hypothesis(estimated.reg,H0)
  # Vergleiche p-value mit Signifikanzniveau
  reject <- ttest[["Pr(>F)"]][2] < sig.level
  
  return(reject)
}

# Funktion, die simulierten t-Test (s.o.) iteriert
# und Anteil der Ablehnungen ausgibt
#
# Argumente:
# ... siehe oben
# Reps    - Iterationen
#
# Funktionswert:
# Anteil der Ablehnungen

rejections.ttest <- function(sig.level=0.05,beta=1,alpha_0=1,alpha_1=0,df_error=100,n=100,Reps=1000)
{
  # Initialisiere reject
  reject <- numeric()
  # Wiederhole Prozedur mit immer neuem seed (i), multipliziere mit 100
  for (i in 1:Reps) reject[i] <- sim.ttest(sig.level,beta,alpha_0,alpha_1,df_error,n,seed=i)*100
  # Bilde Mittelwert (Ablehnungswahrscheinlichkeit...)
  return(mean(reject))
}

#################################################
### (S3)-Klassen (elegant)
#################################################


# S3-plot Methode für Klasse "xfun"
plot.xfun <- function(x)
{
y <- x$f(x$x)
plot(x$x,y,type="l",main=x$name,xlab="x")
}

QuadF <- list(x=seq(-1,1,0.01), f=function(x) 5+0.6*x-0.2*x^2, name="Quadratic Function" )
plot(QuadF)

class(QuadF) <- "xfun"
plot(QuadF)


SinF <- list(x=seq(-1,1,0.01),f=function(x) sin(x), name="Sinus Function" )
class(SinF) <- "xfun"
plot(SinF)




#################################################
### Aufgaben J3: Effizienz
# 1. Lesen Sie Ligges, Kapitel 5

# 2. Gehen Sie Ihre bisherigen Codes (Aufgaben, eigene Codes, ....) durch.
#    Sind diese hinreichend kommentiert?
#    Versuchen Sie diesen zu optimieren: Kommentieren, Profiling, ...