#################################################
### Numerische Optimierung
#################################################

# In VWL, Ökonometrie und Statistik ist häufig Maximierung
# oder Minimierung von Funktionen gefordert. Bei einigen Funktionen 
# (statistische Zielfunktionen wie Likelihood oder Residuenquadratsumme 
# nichtlinearer Regressionen) ist das Optimierungsproblem analytisch 
# nicht lösbar, es existiert keine (algebraische) Formel für den Schätzer (wie etwa bei OLS)

# Ausweg: Es ist möglich, die Funktion NUMERISCH bzw. allgemeiner ITERATIV (in mehreren Schritten)
# zu optimieren. Dazu muss i.A. nur der Funktionswert evaluierbar sein. 

# Achtung: Bei numerischen Algorithmen muss auf Eindeutigkeit der Nullstelle, des Optimums, 
# des Fixpunkts usw geachtet werden. U. a. können unterschiedliche Startwerte zu 
# unterschiedlichen Ergebnissen führen (s. unten)

# Funktionen optimize (univariat), nlm() und optim() (multivariat)
# -> siehe Fortgeschrittene Ökonometrtie für die Darstellung der Methoden
# Hier nur optim verwendet (am allgemeinsten!)
?optim

#### Beispiel 1: Einfache univariate Funktion

f <- function(x) 0.2*(x-1)^4-0.7*x^2-0.3*x+5
x <- seq(0,4,0.01)
plot(x,f(x),type="l")

f(2)
# Minimum bestimmen
# par = Startwert, fn=functionname der Zielfunktion, lower = Untergrenze bei Brent, 
# upper = Obergrenze
?optim

optim(par=NA, fn=f, lower=-10, upper=10, method="Brent")  

abline(v = 2.726503)


### Aufgabe H1
# Modifizieren Sie die Funktion f, so dass vor der Ausgabe des Funktions-
# werts mit der print-Funktion der x-Wert und der Funktionswert ausgegeben
# wird. Führen Sie nun die Optimierung durch. Wie erklären Sie sich das?


### Aufgabe H2 f(x) = exp(x) + log(x)
# Lesen Sie sich die Wikipediaeinleitung zur Bisektion durch. 
# Schreiben Sie nun einen Algorithmus, der Nullstellen einer Funktion auf einem Inter-
# vall [-4,6] sucht und testen Sie ihn an der Funktion 3*x-4
#
# Hinweis:
# Schreiben Sie eine while-Schleife mit der Abbruchbedingung, dass der
# Unterschied zweier Iterationsschritte kleiner 0.01 ist. 
# Berechnen Sie in der Schleife den Mittelpunkt des Intervalls, 
# vergleichen Sie den Wert des Mittelpunkts mit den 
# beiden anderen Werten der Intervallgrenzen, wählen Sie nun für den nächsten Schritt das Intervall, 
# bei dem sich die Nullstelle befinden muss -> Vorzeichenwechsel!



#### Wieso Optimierung in der Ökonometrie?
# Viele Schätzer basieren auf der Minimierung einer Funktion von den Residuen
# Beispiel: Nichtlineare Regression (siehe Fortgeschrittene Ök.)

## Regressionsmodell mit zwei Parametern a1 und a2
a <- c( 1 , 0.6 )
n <- 200

x <- rchisq(n,df=3)
u <-  rnorm(n,sd=0.5)

# Nichtlineares Modell y = a x^b+u
y <-  a[1]*x^a[2] + u

# Meist verwendet man stattdessen log(y) = c0 + c1 log(x) + e
DATA <- data.frame(y,x)
plot(DATA$x,DATA$y)


## Residuenquadratsumme zu minimieren
RSS <- function(a,DATA){
  u_tilde <- DATA$y - a[1]*DATA$x^a[2]
  return(sum(u_tilde^2))
}

## Zielfunktion anschauen für verschiedene Werte von a
length.grid <- 50
a1_grid <- seq(-1,3,length.out=length.grid)
a2_grid <- seq(0,1.5,length.out=length.grid)
a_grid <- as.matrix(expand.grid(a1=a1_grid,a2=a2_grid))

RSS_grid <- apply(a_grid,1,RSS,DATA=DATA)


RSS_grid <- matrix(RSS_grid,length.grid,length.grid)

# Nicht viel zu sehen
contour(a1_grid,a2_grid,RSS_grid,xlab="a1",ylab="a2",nlevels=100)
# -> Logarithmieren
contour(a1_grid,a2_grid,log(RSS_grid),xlab="a1",ylab="a2",nlevels=100)
# 3D-Plot
?persp
persp(RSS_grid,xlab="a1",ylab="a2",phi=10,theta=-40,ticktype="detailed")

# damit das Minimum später auch sinn macht! 
persp(log(RSS_grid),xlab="a1",ylab="a2",phi=10,theta=-40,ticktype="detailed")

## Numerische Minimierung
# par = Startwert
# fn = Zielfunktion

est <- optim(par=c(1,1), fn=RSS, DATA=DATA)   

points(x = est$par[1], y = est$par[2] , col = "red" , pch = 20)
# Fazit: Am Tripel (0.99, 0.61, 45.14) = (a1, a2, f(a1,a2))
# liegt das Minimum der Residuenquadratsumme, da am Punkt (a1, a2)
# die Funktion RSS das Minimum, den geschätzten Wert 45.14 annimmt
# Allgemein:
# Hat man eine Funktion f: R^m -> R^n, so hat der Graph Dimension n+m
# In diesem Fall: R^m=R^2 und R^n=1, also dreidimensional 

?optim

## Bei nichtlinearer Regression einfacher: nls()
est_nls <- nls(y~a1*x^a2, start=c(a1=1,a2=1),data=DATA)
summary(est_nls)




#### Beispiel 3: Startwerte sind wichtig!

# Ziel: Maximieren der Funktion f
f <- function(x) (2*cos(x[1]/0.5)-0.3*x[1]^2)+cos(x[2]^2)-0.5*(x[2]-2)^2

# max statt min: control-Argument benutzen (oder g=-f minimieren)
optim(par=c(1,1),f,control=list(fnscale=-1))
optim(par=c(-0.5,3),f,control=list(fnscale=-1))   # Andere Startwerte...


# Funktion betrachten
l.grid <- 50
grid.x <- as.matrix(expand.grid(x1=seq(-5,5,length.out=l.grid),x2=seq(-5,5,length.out=l.grid)))
grid.f <-apply(grid.x,1,f)
matrix.grid.f <- matrix(grid.f ,l.grid)

image(matrix.grid.f,xlab="x1",ylab="x2")
contour(matrix.grid.f,nlevels=100,add=TRUE)

persp(x = seq(-5,5,length.out=l.grid),
      y = seq(-5,5,length.out=l.grid),
      z = matrix.grid.f,
      xlab = "x1", ylab = "x2", theta = 60, ticktype = "detailed")

# Und Startwerte geschickt wählen
(x.start <- grid.x[which.max(grid.f),])
optim(par=x.start,f,control=list(fnscale=-1))




### Aufgabe H3
# Maximieren Sie die Funktion f
f <- function(x) (2*cos(x[1]/0.5)-0.3*x[1]^2)+cos(x[2]^2)-0.5*(x[2]-2)^2
# Gehen Sie davon aus, dass das Maximum (x1, x2) -5<x1<5 und -5<x2<5 erfüllt.
# Achten Sie auf die Startwerte!




#### Beispiel 3: Maximum Likelihood Schätzung GARCH-Modell (siehe Multivariate
#### Financial Volatility Modelling, Prof. Weber)

# x[t] = h[t]*e[t]
# e[t] ~ IID(0,1)                 'White Noise'
# -> Var(x[t]|I[t-1]) = h[t]^2     Bedingte Varianz ändert sich (ist von Vorperioden abhängig:)
# h[t]^2 = c + g[1]*h[t-1]^2 + ... + g[q]*h[t-q]^2 + a[1]*x[t-1]^2 + ... +  a[p]*x[t-p]^2


### Definiere die Log Likelihood Funktion
### theta = Parametervektor
### x = Vektor (univariate Zeitreihe)

LogLik.GARCH <- function(theta, x)  
{
  c <- theta[1]	# Konstante in GARCH Gleichung
  a <- theta[2]	# ARCH Parameter
  g <- theta[3]	# GARCH Parameter
  # Unzulässige Parameterwerte ausschließen (-> Inf als Wert)
  if (sum(a)+sum(g)>1) return(NA) # Instabil -> nicht zugelassen
  if (any(theta<0)) return(NA)      # Nichtnegativität verletzt -> nicht zugelassen
  
  n <- length(x)			# Stichprobengröße
  h_sq <- rep(var(x),n)	# Vektor der bedingten Varianzen (am Anfang: unbedingte Var.)
  ll <- numeric(n)		# Vektor der (MINUS) Log Likelihoods pro Beobachtung
  
  # Iterativ Bedingte Varianz und minus Likelihood berechnen
  for (i in 2:n)
  {
    h_sq[i] <- c+a*x[i-1]^2+ g*h_sq[i-1]
    ll[i] <- -0.5*(log(2*pi)+log(h_sq[i])+x[i]^2/h_sq[i])     # minus Likelihood
  }
  
  # Funktionswert = Log Likelihood
  return(sum(ll))
}




### Aufgabe H4
# a) Suchen und finden Sie eine möglichst lange Zeitreihe zu täglichen Preisen
#    von Getreidefutures.

# b) Finden Sie für x='Vektor der log-Preisänderungen' (Renditen)
#    den Parametervektor, der die oben definierte Funktion maximiert
#    -> Maximum Likelihood Schätzer

# c) Berechnen Sie die bedingten Varianzen (h_sq)

# d) Plotten Sie untereinander Renditen und bedingte Varianzen
    
# e) Prognostizieren (berechnen) Sie die HEUTIGE bedingte Varianz der Renditen!
# Heute sei Aug 2012, letztes vorhandener Wert von Mar 2012 ->  5-step forecast

# f) Ist diese größer/kleiner als "üblich"?




### Aufgabe H5
# Im File FunktMIN.RData sind die Vektoren y, x und z gespeichert. 
# Betrachten Sie das Regressionsmodell 
# y = beta[1]+beta[2]*x+beta[3]*z^beta[4] + u

                                        
# a) Geben Sie eine Funktion an, die aus den Daten und beliebigen Werten von beta
#    die Residuenquadratsumme berechnet.

# b) Berechnen Sie den Wert von beta, der diese SSR minimiert. 

# Startwerte aus linearen Regression (mit beta[4]=1)

# c) Versuchen Sie, die Funktion f im File FunktMIN.RData
# zu minimieren. Geben Sie auch die Hessematrix an.